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1.0  INTRODUCTION 


The  extraction  of  time-domain  information  from  a  waveform  (e.g.,  pulse  shape) 
can  be  facilitated  by  preprocessing  the  waveform  with  a  “polyline  algorithm,”  i.e., 
an  algorithm  that  approximates  the  waveform  as  a  concatenated  sequence  of  straight- 
line  segments.  The  polyline  approximation,  if  properly  accomplished,  smooths  in¬ 
consequential  noise  structures  while  leaving  significant  structure  well  defined  and 
unsmoothed.  Moreover,  the  polyline  representation  generally  provides  a  large  de¬ 
gree  of  data  compaction  relative  to  the  original  time  series. 

Polyline  approximation  has  recently  been  applied  to  a  great  diversity  of  problem 
areas  (see  the  Bibliography). 

An  early  polyline  application  that  in  many  ways  anticipated  our  present  work  is 
provided  by  McAulay  and  Denlinger1  in  their  development  of  decision-directed 
tracking  algorithms.  An  early  scan-along  algorithm  was  developed  by  Tomek,  based 
on  the  notion  of  finding  the  longest  approximating  line  segment  that  can  be  con¬ 
fined  completely  within  a  pair  of  parallel  lines  separated  by  a  prespecified  error  toler¬ 
ance.2''  Some  of  the  difficulties  experienced  with  Tomek’s  algorithm  are  alleviated 
by  a  scan-along  algorithm  variously  referred  to  as  the  cone  intersection  and  mini¬ 
mum  perimeter  polygon  method,  discovered  independently  by  Williams4  and  by 
Sklansky  and  Gonzalez.' 

Pavlidis1  and  Pavlidis  and  Horowitz6  describe  a  split-and-merge  technique  that 
progressively  improves  on  an  initial  segmentation  until  an  a  priori  error  specifica¬ 
tion  is  satisfied.  More  recently,  Pavlidis7  has  devised  an  algorithm  that  recasts  his 
earlier  split-and-merge  approach  into  a  scan-along  structure  that  he  refers  to  as  a 
hop-along  algorithm. 

The  scan-along  and  split-and-merge  algorithms,  while  effective  in  applications, 
are  not  intended  to  be  optimal  in  any  sense.  Vandewalle*  has  provided  an  algorithm 
that,  while  slow  in  execution,  is  intended  to  provide  an  approximation  that  is  op¬ 
timal,  in  the  sense  of  requiring  the  minimum  number  of  breakpoints  (or  knots)  to 
achieve  a  prespecified  error  norm.  Similar  concepts  are  cast  into  a  somewhat  more 
formal  setting  by  McLaughlin  and  Zacharski,4  who  refer  to  their  algorithm  as  the 
method  of  E-maximal  knots. 

A  dynamic  programming  algorithm  for  polyline  approximation  has  recently  been 
devised  independently  by  Papakonstantinou 10  and  by  Dunham."  These  authors 

'R.  .1.  McAulay  and  E.  J.  Denlinger,  "A  Decision-Directed  Adaptive  Tracker.”  IEEE  Trans.  /I er¬ 
as/?.  Electron.  Syst.  AES-9,  229-236  (1973|. 

:l.  Tomek,  "Two  Algorithms  for  Piecewise  Linear  Continuous  Approximation  of  Functions  of 
One  Variable,”  IEEE  Trans.  Coinpitt.  0-23,  445-448  (1974). 

'T.  Pavlidis,  Structural  Pattern  Recognition,  Springer-Verlag,  Berlin  (1977). 

JC.  M.  Williams,  “An  Efficient  Algorithm  for  the  Piecewise  Linear  Approximation  of  Planar 
Curves,”  Coinput.  Graph.  Image  Process.  8,  286-293  (1978). 

'.I.  Sklansky  and  V.  Gonzalez,  “Fast  Polygonal  Approximation  of  Digitized  Curves,”  in  Proc. 
1979  IEEE  Computer  Society  Conf.  Pattern  Recognition  Image  Processing,  pp,  604-609  (1979). 
"T.  Pavlidis  and  S.  I  .  Horowitz,  "Segmentation  of  Plane  Curves,"  IEEE  Trans.  C'oinpul.  0-2.3, 
860-870  (1974). 

T.  Pavlidis,  Algorithms  for  Graphics  and  Image  Processing.  Computer  Science  Press,  Rockville, 
MJ.  (1982). 

\l.  Vandewallc,  “On  the  Calculation  of  the  Piecewise  I  inear  Approximation  to  a  Discrete  Func¬ 
tion,"  IEEE  Trans.  C'oinpul.  0-24,  843-846  (1975). 

S’H.  W.  McLaughlin  and  .(.  J.  /acharski,  “Segmented  Approximation,"  in  Approximation  The¬ 
ory.  E.  W.  Cheney,  cd..  Academic  Press.  New  York.  pp.  647-654  (1980). 

"'G.  Papakonstantinou.  “Optimal  Polygonal  Approximation  of  Digital  Curves,”  Signal  Process. 
8,  131-135  (1985). 

"J.  G.  Dunham,  “Optimum  Uniform  Piecewise  Linear  .Approximation  of  Planar  Curves,"  IEEE 
frans.  Pattern  Anal.  Mach.  Intel!.  PAMI-8.  67-75  (1986). 
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claim  optimaliiv  tor  their  algorithms,  in  the  sense  of  obtaining  an  approximation 
with  a  globally  minimum  number  of  knots,  for  a  prespecified  peak  fitting  error. 
The  dynamic  programming  algorithms'"  1  call  upon  the  cone  intersection  method 
as  a  subroutine. 

Numerical  experience  with  the  "optimal”  algorithms  indicates  that  the  dynamic 
programming  algorithm  prosides  both  reduced  fitting  error  and  faster  execution 
speed,  relative  to  the  method  of  L-maximal  knots,  as  discussed  in  Section  6  below. 

In  addition  to  providing  a  new  dynamic  programming  algorithm,  Dunham  presents 
performance  comparisons  of  sev  eral  prior  polyline  algorithms  against  three  test  con¬ 
tours."  The  measures  of  algorithm  performance  used  by  Dunham  are  execution 
speed,  degree  of  data  compaction,  and  peak  absolute  fitting  error;  algorithms  evalu¬ 
ated  included  those  by  Pavlidis  and  Horowitz,"  I’avlidis.  Williams,4  Sklanskv  and 
Gonzalez,'  Uadi ’ i  and  Peikari."  Rainer,''  Roberge. and  Dunham."  Dunham’s 
results  show  that: 

1.  Roberge’s  algorithm,  while  consistently  fastest  in  execution,  typically  achieved 
only  half  the  data  compaction  ratio  of  the  optimal  dynamic  programming  so¬ 
lution  (for  fixed  peak  error): 

2.  Tlte  Sklanskv.  Gonzalez  algorithm  executed  about  as  fast  as  the  fastest  of  the 
other  algorithms  (not  including  Roberge’s)  and  consistently  achieved  data  com¬ 
paction  ratios  almost  as  good  as  the  dynamic  programming  algorithm. 

Consequently,  the  Sklanskv  ''Gonzalez,  algorithm  may  be  considered  a  benchmark 
against  which  to  compare  the  performance  of  other  methods. 

Pavlidis,  in  editorial  commentary,  has  noted  that  performance  comparisons 
such  as  Dunham's  are  clouded  bv  uncertainties  in  algorithm  implementation.  How¬ 
ever,  also  as  noted  by  Pavlidis.  such  numerical  studies  and  algorithm  intercompari¬ 
sons  arc  nonetheless  valuable,  following  Dunham,  we  presently  both  introduce  new 
algorithms  and  compare  them  to  prior  ones.  (  omparing  our  work  with  Dunham’s: 

1.  Our  interest  is  in  approximating  waveforms  i  at  her  than  two-dimensional 
contours. 

2.  We  have  expanded  the  scope  of  etror  metrics  to  include  average  error  (i.e., 
bias),  root  mean  square  error,  and  average  absolute  error. 

3.  Our  selection  of  polvlinc  algorithms  tot  evaluation  and  comparison,  partially 
overtopping  Dunham's,  includes  those  of  Pavlidis,  Papakonstantinou, 111 
Sklanskv  and  Gonzalez.  Wall  and  Danielsson, "  Tomek,"  and  McLaughlin 
and  Zacharski. '  as  well  as  the  new  algotithms  presented  in  this  report. 

As  discussed  in  Section  6,  the  relative  performance  of  the  various  algorithms  turns 
out  to  be  somewhat  dependent  on  the  error  metric.  In  particular,  the  original  al¬ 
gorithms  presented  hete  are  considerable  superior  to  prior  algorithms,  with  respect 
to  bias  and  room  mean  square  etror. 

Our  new  algorithms  are  of  two  varieties,  which  we  refer  to  as  recursive  least-squares 
(RI.S)  and  generalized  likelihood  ratio  |GI  R)  algorithms. 


T.  Badi'i  ami  It.  Petkait.  "I  nncuuit.il  Xppiu'.iip.ation  el  Planar  C  lines  via  Adaptive  Segmenta¬ 
tion."  Ini.  ./.  Sts/en/s  Si  i.  13.  Mr  fi'4  1 1 '*s2 ) . 

'U.  Ranter,  "An  itetative  Ptoceditii  I  n  ih,  I'olie.iiial  \ppro\imatiim  ol  Plane  C  urves."  Ctnn- 
pul.  Graph.  Imam'  /Voces'.  I.  244  2~4»  (l'P2t 

iJ.  Roberge.  "A  Data  Reduction  Xleortihtn  lot  Planar Wanes. “  (inn/uti  \  i\inn  Graph.  Image 
Process.  2*/.  UsN  Ids  tlUSS). 

I.  Pavlidis.  "Pdiliirt.il-  Papcn  on  shape  Snahiis,"  //  Z  /  him  s.  Pal  lari)  \nal.  Mach.  Ink'll. 
PAVII-fi.  I  (PPr.), 

"Is.  Wall  and  P.  l)anieh»on.  "  \  I  a-!  Sequential  Method  lot  Polvgonal  Approvintation  of  Digi¬ 
tized  Curses."  f  ampul  Graph.  Imam  Prints.  2X.  22u  22~  1 1 ')S4 ) 
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Our  RLS  algorithms  are  derived  in  Section  4  based  on  a  simple  Kalman  filter. 
Our  numerical  experience  with  these  algorithms  is  that,  relative  to  the 
Sklansky/Gonzalez.  algorithm,  R1  S  has 

•  Faster  execution  speed 

•  Slightly  worse  peak  error 

•  Less  average  error 

•  Smaller  room  mean  square  error. 

The  original  application  of  the  statistical  GL.R  method  to  formulating  the  prob¬ 
lem  of  change  detection  in  linear  systems  was  provided  by  McAuIay  and  Denlinger. 1 
A  general  solution  to  the  GL.R  formulation  was  provided  by  Willsky  and  Jones.1 
In  Section  2  of  this  report  we  review  the  Willskv/Jones  GLR  formalism  that,  in 
Section  5,  we  particularize  to  the  problem  of  polyline  approximation. 

The  most  direct  antecedent  of  our  work  is  McAuIay  and  Denlinger. 1  Features 
common  to  our  work  and  theirs  are: 

1.  Use  of  a  GLR  formulation  for  change  detection; 

2.  “...A  piecewise  linear  model...  to  (describe)  the  noiseless  evolution...”  of  mea¬ 
sured  waveforms;1 

3.  Tracking  of  “...straight-line  segments  with  a  simple  two-term  Kalman 
filter...”. 1 

Some  points  of  distinction  between  our  approach  and  theirs  are: 

1.  We  develop  a  one-parameter  recursive  regression  as  our  candidate  no-jump 
Kalman  filter.  The  corresponding  filter  used  in  Ret.  1  is  unknown  since  its 
discussion  is  quite  brief. 

2.  Our  “jump  signature”  solution  in  Section  5.3  is  exact,  derived  as  a  special 
case  of  the  general  solution  in  Ref.  17.  The  corresponding  signature  presented 
in  Ref.  1  is  approximate,  being  based  on  numerical  experience  (cf.  the  discus¬ 
sion  in  connection  with  Eq.  106). 

Although  Ref.  I  is  widely  cited  as  the  original  GLR  change  detection  formula¬ 
tion,  the  casting  of  target  tracking  problems  in  terms  of  an  equivalent  polyline  ap¬ 
proximation  has  apparently  gone  unnoticed  in  the  literature,  perhaps  due  to  the 
brevity  of  this  part  of  the  discussion  in  Ref.  1.  Nevert holes  ,  we  see  this  as  poten¬ 
tially  of  high  interest,  since  algorithms  for  polyline  approximation  can  be  applied 
to  target  tracking  and  vice  versa.  While  there  has  very  recently  been  some  renewed 
interest  in  applying  polyline  approximation  methodology  to  problems  in  radar  track¬ 
ing,1'’  there  may  still  be  considerable  unexploited  potential  for  cross-fertilization  be¬ 
tween  the  two  areas. 


1  A.  S.  Willsky  and  H.  I  .  Jones.  “A  Generalized  I  ikelihood  Ratio  Approach  to  the  Detection 
and  Estimation  of  Jumps  in  l  inear  Systems,”  //  /  /:  Trans.  Amo.  Control  AC-21,  108-1 12  (1976), 
“S.  I  Haase.  “Advanced  Radar  Tracking  Techniques,”  in  IR&D/H&P  Program  Plan  Vol.  II 
PYI9SS.  IHC  API  ,  1  aurcl.  Mil.,  pp.  li-J'S  to  11-397  (1988). 
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2.0  DETECTING  CHANGES  IN  LINEAR  SYSTEMS 


2.1  PROBLEM  FORMULATION 

This  section  follows  Willsky  and  Jones'7  and  like  them  provides  results  without 
derivations.  For  completeness,  the  appropriate  derivations  are  provided  in  Appen¬ 
dix  A. 

We  assume  a  discrete-time  system  whose  dynamics  evolve  according  to  the  state 
and  measurement  equations: 

.v(A-  +  1 )  =  sJ5 ( At  +  1  ,k)  ■  x(k)  4-  T(A)  •  w(k)  +  ■  v  ,  (1) 

z(k  +  1)  =  H(k  +  1)  ■  ,v(A-  +  1)  +  v(A-  +  1)  ,  (2) 


with  all  quantities  being  real  valued  (see  the  list  below  for  definitions).  Quantities 
v,  w,  i>,  z,  and  v  may  be  vectors,  while  4>,  P,  and  H  are  real  matrices  of  commen¬ 
surate  dimensionality. 

Quantity  Definition 

k  Discrete  time 

.v  State  vector 

if  State  noise 

c  Measurement /observation 

v  Measurement  noise 

6  Dirac  delta  function 

i>  Jump  amplitude 

0  Jump  occurrence  time 

4>  State  transition  matrix 

H  Measurement  matrix 


The  noise  sequences  \v(k)  and  v(A)  are  assumed  to  be  zero-mean  and  Gaussian, 
with  covariances 

cov| >v(A'))  =  £[w(A-)w'(A-)|  =  Q(k)  .  (3) 

cov[v(A))  =  R(k)  .  (4) 

where  in  Eq.  3  £(  ■]  denotes  ensemble  expectation,  and  a  primed  quantity  denotes 
the  vector  transpose. 

The  Dirac  delta  term  in  Eq.  1  indicates  the  presence  of  a  discontinuity  (or  “jump”) 
in  what  would  otherwise  be  a  smooth  temporal  evolution  in  the  state  vector.  The 
jump  amplitude  r  and  time  of  occurrence  6  are  presumed  to  be  unknown,  subsum¬ 
ing  the  possibility  that  a  jump  never  occurs  (i>  =  0  or  0  -  oo ). 

The  objective  of  jump  detection  processing  is  to  determine  whether  a  jump  has 
occurred  and,  if  so,  to  establish  accurate  estimates  for  v  and  0  by  means  of  process¬ 
ing  operations  performed  on  the  measurements  '(A). 

Since  Eq.  1  has  the  stochastic  driving  term  vv,  the  state  x  must  necessarily  be  a 
random  process. 
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2.2  KALMAN  PREPROCESSOR 

Our  approach  to  detecting  the  presence  of  a  jump  in  E;,.  1  is  to  process  the  data 
with  a  filter  that  will  be  optimal  if  no  jump  occurs  and  to  inoniior  the  filter’s  per¬ 
formance  to  assess  the  continuing  validity  of  the  no-jump  hypothesis.  In  the  control 
systems  literature,  this  approach  is  sometimes  referred  to  as  analytical  redundancy, 
and  the  filter  is  sometimes  called  a  no-fail  observer  or  normal-mode  observer. 19 

Subject  to  the  no-jump  hypothesis  (i.e.,  prior  to  the  jump,  when  k  <  8),  the  op¬ 
timal  approach  for  estimating  the  state  is  provided  by  a  Kalman-Bucy  filter: 

A-(0|0),  initial  state  estimate  (5) 

P(0|0),  initial  covariance  (6) 

x(k  +  1 1 A)  =  4>(Ar  +  l.Ar)  •  ,v(A|A)  (7) 

z(k  +  1 1 A:)  =  H(k  +  1 )  •  x(k  +  1 1 A)  (8) 

7 (k  +  1)  =  z(k  +  1)  -  z(k  +  1 1 A)  (9) 

P(k  +  1 1 A)  =  4>(A  +  1  ,k)  ■  P(k\k)  ■  4 >'<A  +  1  ,k) 

+  T(A')  •  Q(k)  ■  r'  (A)  (10) 

V(k  +  1)  =  H(k  +  1)  •  P(k  +  1 1 A)  •  H'(k  +  1) 

+  R(k  +  1)  (11) 

K(k  +  1)  =  P(k  4-  1 1 A)  ■  H'  (k  +  1)  •  V~'(k  +  1)  (12) 

.v( A  +1[  +  1)  =  x(k  +  1 1 A)  +  K(k  +  1)  •  7(A  +  1)  (13) 

P(k  +  1 1 A  +  1)  =  [/  -  A'(A  +  1)  •  H(k  +  1)]  ■  P(k  +  1 1 A)  (14) 

with  all  quantities  in  Eqs.  5-14  being  real  valued.  (See  the  list  below  for  definitions; 

note  that  quantity  Z1  denotes  the  sequence  of  observations  c(l),  z(2) . z(j).)  The 

flow  of  operations  in  the  Kalman  filter  is  depicted  in  I  ig.  1. 


Quantity 

Definition 

-v  ( A 1 7 ) 

Conditional  mean  of  x(k),  given  Z' 

z(k\j) 

Conditional  mean  of  z{k),  given  Z' 

7(A) 

Innovation 

P(k\j) 

Covariance  of  x(k\j) 

Vik) 

Covariance  of  7(A) 

K(k) 

Kalman  gain 

’’’A.  Madiwale  and  B.  I  ricdland,  “C'oinaprison  of  Innovations-Rased  Analytical  Redundancy 
Methods.”  in  Prnc.  198.1  Am.  Control  Con/.,  pp.  940-945  (19X3). 
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Figure  1  Kalman  filter  flow  of  operations. 

2.3  ca  R  KOKMAI  ISM 

As  first  noted  by  McAulay  and  Denlinger.1  the  linearity  of  the  system  equations 
allows  their  solution  in  response  to  a  jump  excitation  to  be  written  as  the  sum  of 
a  jump-independent  stochastic  component  (.v, )  and  a  second  component  (ax)  linear¬ 
is  proportional  to  the  jump  amplitude.  Moreover,  with  the  jump  amplitude  and 
time  taken  as  deterministic  (though  unknown)  quantities,  the  jump-dependent  com¬ 
ponent  ,v:  is  seen  to  be  deterministic. 

Following  Ref.  17,  we  write  the  state  as 

v ( A )  -  .v,  (A)  +  4>(A,ff)  •  e  ,  (15) 

where  x,  (A)  is  the  value  that  v(A)  would  take  in  the  absence  of  a  jump  (i<  =  0), 
and  4>(A ,0)  •  r  is  the  state  perturbation  in  response  to  a  jump,  where 

<Hk.li)  -  0,  k  <  D  (16) 

<YWJ>)  I  (17) 

<1>(A  *  IJf)  -  <J>(A  t  I. A)  ■  <J>(A,0)  .  (18) 

I  he  proof  of  Fqs.  16-18  is  prov  ided  in  Appendix  A.  along  with  the  proofs  for  oth¬ 
er  results  presented  in  this  section. 

From  Fqs.  2  and  15. 

-MA  I  =  (A)  +  //(A)  ■  <I>(A,W)  •  f -  .  (19) 

Similarly,  the  linearity  of  the  Kalman  filter  equations  allows  their  solution  to  be 

written  in  a  form  analogous  to  Fqs.  15  and  19,  vi/.,  for  the  state  estimate, 

v< A  A)  =  a,  (A: A)  +  /  ( kj> )  ■  r  ,  (20) 
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and  for  the  innovation. 


y(k)  =  7l  (*)  +  G(k-fi)  ■  »  . 


The  quantities  F'(k;6)  and  G(k\Q )  appearing  in  Eqs.  20  and  21  are  sometimes  referred 
to  as  the  jump  “signatures”;  they  characterize  the  Kalman  filter’s  response  to  a  jump 
excitation  in  the  system  equations. 

It  clearly  must  be  true  that 

F(k-fi)  =  0  0 

G(k;8)  =  0  J  A  <  d  (22) 

in  order  to  preserve  the  definitions  of  .v,  and  7,  as  the  Kalman  filter  responses  in 
the  absence  of  a  jump. 

It  can  be  shown  that  G(k;6)  may  be  written  in  terms  of  F(k;6): 

G(k;8)  =  H(k)  ■  [$(*■, 0)  -  $(k,k  -  1)  •  F(k  -  1;0)]  .  (23) 

The  function  F(k\8)  appearing  in  Eqs.  20  and  23  is  given  by 


F(k;6)  =  £  Q(k-J)  ■  K(J)  •  H(J)  ■  4>O,0)  ,  (24) 

j=» 

where  the  auxiliary  variable  Q(k  J)  is  obtained  recursively  as 

9(k;0)  =  [I  -  K(k)  ■  H(k )]  ■  *(*,*  -  1)  ■  Q{k  -  1,0)  (25) 

9(k;6)  =  0,  k  <  6  (26) 

€3(0,0)  =  /  (27) 

Assuming  that  the  measurements  z{k)  generated  according  to  Eqs.  1  and  2  are 
processed  by  the  Kalman  filter,  Eqs.  5-14,  the  innovation  y{k)  will  be  a  zero-mean 
Gaussian  process  (7, )  until  the  jump  occurs  at  k  =  6,  whereupon  y(k)  will  develop 
a  bias,  G(k\Q)  ■  v.  Jump  detection  is  accomplished  by  detecting  the  innovation  bias, 
which,  in  turn,  is  accomplished  by  a  form  of  matched  filtering  applied  to  the  inno¬ 
vation  process. 

Writing  the  GLR  estimates  for  v  and  0  as  P  and  0,  it  can  be  shown  that 

Hk)  =  C  '(k;9)  ■  d(k;6)\„ =(j(k)  ,  (28) 

where  d{k\6)  is  obtained  via  a  matched  linear  filtering  operation  on  the  innovation 
process, 


d(k\0)  =  £  G’(n\8)  ■  V  '  (n)  ■  7  (n)  , 


and  matrix  C  '(k;8)  is  the  covariance  of  P.  computed  as 


C{k:8)  =  £  G’  (n;0)  •  V  '  (n)  ■  G(rr,8) 
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The  GLR  estimate  for  jump  time,  6 ,  is  obtained  as  the  value  of  d  that  maximizes 
the  generalized  likelihood  ratio,  i.e., 

I \k;0)  =  max  t(k\d)  ,  (31) 

u  i.: . * 


Kk;d)  =  d'(k\6)  ■  C  \k;8)  ■  d(k;6)  .  (32) 

The  GLR  algorithm  determines  that  a  jump  has  occurred  whenever  t(k;d)  exceeds 
a  fixed  threshold,  e;  the  GLR  estimate  for  the  jump  occurrence  time  is  then  provid¬ 
ed  by  Eq.  31,  and  the  GLR  estimate  for  jump  amplitude  v  is  given  explicitly  by  Eq.  28. 

We  note  that  the  function  C(k\d)  ghen  by  Eq.  30  is  computable  off  line,  as  is 
the  function  [G'(/';0)  •  V  '(/)]  in  Eq.  29. 

The  GLR  filter  flow  of  operations  is  diagrammed  in  Fig.  2a.  Maximization  in 
Eq.  31  occurs  over  k  values  of  0,  corresponding  to  the  k  parallel  branches  in  Fig. 
2a.  Because  the  number  of  calculations  required  to  implement  Eq.  31  increases  with 
time,  implementing  “full  GLR”  requires  progressively  increasing  numbers  of  cal¬ 
culations  at  successive  time  steps — an  undesirable  property.  An  approach  to  bounding 
the  number  of  required  calculations  is  to  restrict  the  range  of  0  values  in  Eq.  31 
to  the  M  “most  recent”  values,  i.e.,  to  the  range  [(k  -  M  +  1  ),Ar]  (Fig.  2b). 

When  the  jump  vector  v  can  be  written  as 


where  a  is  an  unknown  scalar  and _/(<))  is  a  known  vector  function  of  0,  the  formula¬ 
tions  for  v  and  t(k\8)  may  be  written  in  the  following  forms,  alternatives  to  Eqs. 
28  and  32: 

v(k)  =  a(k;8)  ■  f(6)  (34) 

at(k;8)  =  b(k;6)/a(k;8)  (35) 

d k-6 )  -  b2(k\0)/a{k\0)  (36) 


a{k\9)  =  /'«?)  •  C(k;8)  ■  f(6) 


b(k\6)  =  f'(6)  ■  d(k\8)  . 


We  note  that  Basseville  and  Benveniste:"  have  suggested  a  modified  form  of  the 
Willsky/Jones  algorithm  in  which  thresholding  is  performed  on  the  estimated  jump 
amplitude  v,  given  by  Eq.  28.  Much  simpler  still,  in  Section  4.2  we  propose  threshold¬ 
ing  the  normalized  innovation  y(k)/( V(k))  directly,  i.e.,  using  the  Kalman  filter 
alone,  with  no  additional  processing  other  than  thresholding,  for  jump  detection. 

:"M.  Basseville  and  A.  Benvenistc,  "Design  and  Comparative  Study  of  Some  Sequential  Jump  De¬ 
tection  Algorithms  for  Digital  Signals,"  IEEE  Trans.  Acoust.  Speech  Signal  Process.  ASSP-31, 
521-535  (1983). 
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Measurements 

z[k) 


(a)  Full  GLR. 


Measurements 

zik) 


(b)  Windowed  GLR  (M  =  20). 


Figure  2  GLR  filter  flow  of  operations. 


The  GLR  algorithm  is  not  altogether  simple  to  understand.  As  perhaps  the  sim¬ 
plest  application  yet  put  forward  of  the  GLR  method,  our  GLR  polyline  develop¬ 
ment  in  Section  5  helps  cast  insight  into  the  structure  of  GLR  generally.  Some  of 
the  simplifying  features  of  our  development  are: 

1.  Low  dimensionality:  2  x  1  vectors  and  2x2  matrices  are  our  highest- 
dimensional  entities. 

2.  Closed  forms:  A  number  of  quantities,  generally  available  only  as  the  solution 
of  recursion  relations,  arc  obtained  in  closed  form  (see  Table  1). 
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Table  1 

Application  of  GLR  to  polyline  approximation  results  in  the  development  of 
closed-form  expressions  for  these  quantities. 


Variable 

Definition 

( )pen  Form 

Closed  Form 

AT(A') 

Kalman  gain 

Lq.  12 

Eq.  62 

HA) 

Innovation  covariance 

Lq.  1! 

Lq.  63 

O(A;0) 

Auxiliary  variable 

Lqs.  25-27 

Lqs.  87  and  89 

/■(A;ff) 

State  signature 

Lq.  24 

Lq.  95 

G(A;«) 

Innovation  signature 

Lq.  23 

Eq.  97 

C(A;fl) 

f  inverse  covariance 

Lq.  30 

Eq.  109 

3.0  STATE  VARIABLE  REPRESENTATION  OF  NOISE-FREE 

POLYLINES 


In  this  section  we  show  that  the  polyline  waveform  approximation  problem  can 
be  given  a  state  variable  formulation,  as  required  for  application  of  the  GLR  for¬ 
malism. 

Underlying  our  development  is  the  assumption  that  waveforms  of  interest  may 
be  approximated  as  noise-free  polylines,  with  additive,  zero-mean,  stationary  Gaussian 
noise.  Thus,  our  interest  in  this  section  is  to  cast  the  equations  of  a  noise-free  poly¬ 
line  into  the  form  of  Lqs.  I  and  2,  in  which  the  noise  terms  have  been  set  to  zero: 


.v(A'  +  1 )  =  (F(A  +  I, A  )  •  ,v(A)  -t-  .  i  •  r;  , 

-(A  +  1)  =  H( k  +  1)  •  ,v(A  +  1)  . 


We  assume  that  the  ordinate  of  the  initial  breakpoint  is  known;  subtracting  the 
initial  ordinate  then  results  in  a  polyline  whose  initial  breakpoint  is  at  coordinates 
(A„,0).  For  further  discussion,  see  Fig.  3,  where  we  define 


/  =  k  A„  , 

(40) 

r  -  0  A , ,  . 

Consideration  of  Fig.  3  shows  that  we  ma>  write  for  the  polyline  slope 

\(i)  n„  •  o  -  /((/  r)  .  (41) 

where  /<(-)  is  the  unit-step  function.  For  offset  we  may  write 

f(/)  K  ■  r  ■  /((/  r)  .  (42) 
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Figure  3  (a)  Noise-free,  two-segment  polyline;  (b)  slope;  (c)  offset. 
Also,  we  note  that 


n(j  ~  r  +  i)  -  n(j  -  r)  =  8JtU  ,  (43) 

where  5  is  the  Dirac  delta. 

From  the  foregoing,  we  obtain  the  desired  state-variable  representation  for  the 
noise-free  polyline: 


(44) 


and 


P 


+  j)  =  i  i  j  i  ■  [$;} 


Comparing  Eqs.  44  and  45  with  Eqs.  38  and  39,  we  obtain 


x(kt)  +  j) 


(45) 


(46) 


THE  JOHNS  HOPKINS  UNIVERSITY 

APPLIED  PHYSICS  LABORATORY 

LAUREL,  MARYLAND 


4>{A'  +  1  ,k)  =  / 


lr\ 


H  ( ku  +  j)  =  |  1  j 


(47) 

(48) 

(49) 


Comparing  Eqs.  48  and  33,  we  see  that  they  are  of  the  same  form,  permitting 
us  to  write 


+  r) 


rv„'H7) 


(50) 


4.0  RECURSIVE  LEAST-SQUARES  ALGORITHMS 


This  section  comprises  three  subsections.  In  4. 1  we  derive  a  recursive  least-squares 
regression  formula  for  fitting  a  straight  line  to  a  set  of  random  data,  as  a  special 
case  of  the  Kalman  filter.  We  assume  that  the  initial  coordinates  of  the  fitted  line 
are  known  and  only  the  line  slope  must  be  estimated.  Since  only  a  single  parameter 
(slope)  is  being  estimated,  our  resulting  “one-parameter  recursive  regression”  is  some¬ 
what  different  from  the  usual  two-parameter  intercept /slope  regression  formula. 21 
In  4.2  we  use  our  one-parameter  regression  formula  to  derive  two  algorithms  for 
polyline  approximation,  which  we  subsequently  refer  to  as  RLS  algorithms.  The  per¬ 
formance  of  these  algorithms  is  compared  in  Section  6  to  other  algorithms.  In  4.3 
we  slightly  recast  our  one-parameter  regression  results  into  a  form  suitable  for  use 
with  the  general  GLR  formulation  presented  in  Section  2. 


4.1  ONE-PARAMETER  RECURSIVE  REGRESSION 

Takings  <  r  in  Eqs.  41  and  42,  we  obtain  the  equations  of  a  straight  line,  viz.. 


.»•(/)  -  0 

(51a) 

v(/)  o„ 

(51b) 

According  to  Eqs.  51,  the  line  passes  through  initial  coordinates  (A„,0),  with  slope 
0(1  • 

From  Eqs.  45  and  51a,  and  adding  a  noise  term  (as  in  Eq.  2),  we  obtain 


:'Y.  Bar-Shalom  and  T.  1  .  I  orimann.  Fracking  and  Dam  Association ,  Academic  Press  (1988). 
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3 


I'M 


z(k)  —  {  ®  >  J  ^  ®  (52) 

Z(  )  "  [7  •  50)  +  v(Ar)  ,  0  <  j  =  (*  -  *„)  <  r  ’ 

where  we  assume  that  v(k)  is  zero  mean,  Gaussian,  and  stationary: 

R(k)  =  cov[v(A:)]  =  a 2  .  (53) 

Under  the  assumption  that  the  initial  ordinate  of  the  line  is  known  and  can  be  sub¬ 
tracted  to  obtain  a  line  that  goes  through  (Ar0 ,0)  (cf.  Fig.  3),  the  single  unknown 
parameter  remaining  is  simply  the  line’s  slope,  s(J)  =  a,,.  Although  our  constraint, 
z(tt0)  =  0,  may  seem  artificial,  it  is  justified  later  in  connection  with  our  applica¬ 
tion  of  the  results  in  this  section. 

A  simple  recursive  procedure  for  estimating  the  parameter  s(j)  =  a0  in  Eq.  52 
may  be  developed  from  the  Kalman  filter  Eqs.  5-14,  in  which  all  generally  vector 
and  matrix  quantities  reduce  presently  to  scalars.  We  identify  the  line  slope  s(J )  as 
the  “system  state,”  which  is  time-independent  (Eq.  51b).  Consequently, 

s(n\j)  =  s(m\J)  ,  (54) 

Ps(n\j)  =  PJm\j)  ,  (55) 

for  all  values  of  n,  m,  and  j,  and  where  we  write  P,  as  the  covariance  of  s.  It  fol¬ 
lows  that  we  may  simplify  notation  by  writing 

s(j)  =  5(w|y)  , 

(56) 

PSU)  *  PAn\j)  . 

Equation  51b  may  be  written  in  the  form  of  Eq.  1  by  making  the  identification 

m  +  l,k)  =  1  (57) 

w(k)  =  0  .  (58) 

From  Eqs.  58  and  3,  Q{k)  =  0.  Similarly,  Eq.  52  can  be  cast  in  the  form  of  Eq. 
2  by  writing 

H(k)  =  (K  —  ka)  ~  j  ■  (59) 

Using  Eqs.  53  and  56-59,  we  obtain,  corresponding  to  the  Kalman  Filter  Eqs.  8-14, 
m  -  i)  =  j  ■  s<j  -  i)  (8') 

y(k)  =  z(k )  —  z(Ar  —  I)  (9') 

PsU)  =  PsU)  U0') 

V(k)  =  j2  ■  P,{j  -  1)  +  a2  (ID 

KSU)  =  j  ■  PsU  ~  \)/V(k)  (12') 

s(j)  =  S(j  -  1)  +  A',0)  •  y(k)  (13') 

psu)  =  [i-;-  k,v))  -  PsU  -  i) .  04  ) 


P 


STV1 


v.  -r, 


WWW 


vv-w-v 
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where  it  is  understood  that  A  =  (A„  +  j)  in  Eqs.  8'  -14'. 

As  shown  in  Appendix  B,  the  desired  recursive  estimation  procedure  for  line  slope 
is  obtained  from  Eqs.  8 '  - 1 4 '  as 

s{J)  =  *0  -  I)  +  A,  0)  ■  7(A)  .  j  =  (k  -  A„)  >  2  ,  (60) 

where 

7(A)  =  z(k)  -  j  ■  s(j  -  1)  ,  (61) 


and 


a;  (./) 


6 

U  +~  1)  ■  ( 2/  +  1) 


(62) 


We  also  find  in  Appendix  B  that 


and 


'(A)  ■  Kd i) 


•M!* 


(63) 


Ps  U) 


6a" 

j  •  0  +  1 )  •  ( 2j  +  1 ) 


(64) 


It  can  be  shown  that  Eqs.  60-64  are  equivalent  to  a  recursive  least-squares  regres¬ 
sion  for  line  slope,  derived  on  the  assumption  that  the  initial  ordinate  is  specified. 
The  recursion  is  initialized  by  taking 

*0)  =  z(k0  +  1)  -  c(A„ )  .  (65) 

The  first  loop  of  recursive  continuation  follows  from  Eqs.  65,  60,  and  61,  in  the 
following  sequence: 

7<A„  +  2)  =  c(A„  +  2)  -  2  •  .9(1)  ,  (66) 

.9(2)  =i(l)  +  A,  ( 2 )  •  7 (A, i  +  2)  .  (67) 

The  slope  estimate  .9(1)  and  measurement  z(k„  +  2)  are  used  in  Eqs.  66  and  67  to 
obtain  the  next  slope  estimate,  .9(2).  More  generally,  we  interpret  Eqs.  60  and  61 
as  a  prescription  for  taking  a  slope  estimate  .v (/  -  I )  and  measurement  c(A)  to  de¬ 
rive  the  next  slope  estimate  s(j). 

Inspection  of  Eq.  62  reveals  two  properties  of  the  Kalman  gain  that  are  generally 
true  even  for  more  complex,  higher-dimensional  Kalman  filters,  viz., 

1 .  The  Kalman  gain,  though  time-dependent,  is  independent  of  the  data  and  may 
be  computed  off  line. 

2.  The  Kalman  gain  becomes  progessively  smaller  as  increasing  amounts  of  data 
are  processed,  i.e.. 
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It  follows  as  a  consequence  of  Eqs.  60  and  68  that  the  RLS  slope  estimate  s(J)  be¬ 
comes  progressively  less  sensitive  to  the  measurements  as  increasing  amounts  of  data 
are  processed.  In  a  sense,  the  filter  becomes  increasingly  “satisfied”  with  the  good¬ 
ness  of  its  slope  estimate.  This  may  be  seen  also  from  Eq.  64:  the  variance  of  the 
estimate,  as  calculated  by  the  filter,  decreases  to  zero.  Naturally,  the  filter’s  assess¬ 
ment  of  its  own  performance,  as  provided  by  Eq.  64,  is  only  valid  if  the  underlying 
waveform  model  is  correct,  i.e.,  if  the  data  realize  a  process  composed  of  a  straight 
line  and  additive,  white,  Gaussian  noise. 

Our  reasons  for  assuming  that  z(k0)  is  known,  and  our  basis  for  choosing  a  val¬ 
ue  for  z(k0),  are  discussed  in  the  next  section. 


4.2  RLS  POLYLINE  ALGORITHMS 

In  this  section  we  describe  two  polyline  approximation  algorithms  based  on  the 
premise  that  the  development  of  a  nonzero  trend  in  y(k)  indicates  that  the  underly¬ 
ing  model,  Eq.  52,  is  no  longer  valid  and  that  the  measurements  can  no  longer  be 
fitted  adequately  with  a  single  straight  line. 

Both  of  our  simple  RLS  polyline  algorithms  are  structured  as  follows  (Fig.  4). 

1.  Equations  60  and  61  are  applied  to  the  data  to  generate  the  innovation  se¬ 
quence  y(k). 

2.  The  innovation  sequence  is  tested  in  some  fashion  for  the  development  of  non¬ 
zero  bias. 

3.  When  a  nonzero  innovation  bias  is  detected,  say,  at  time  k  =  (8  +  1)  = 
(k0  +  r  +  1),  a  breakpoint  or  “knot”  is  introduced  in  the  waveform  approx¬ 
imation  at  time  k  -  0.  The  ordinate  of  the  knot  is  estimated  as 

Z(0)  =  (r  +  1 )  ■  i(r)  ,  (69) 

which  is  obtained  by  substituting  (k  -  1)  =  d  =  (kn  +  r )  into  Eq.  8',  and 
s(r)  is  the  last  slope  estimate  generated  by  the  RLS  algorithm. 


Figure  4  Structure  of  RLS  polyline  algorithms. 
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4.  The  RLS  algorithm  is  reinitialized  hv  making  the  assignments 
z(k)  -  ;(A)  -  zV»  .  k  >  0 


The  purpose  of  Eq.  70  is  to  translate  the  polyline  such  that  the  origin  of  the  new 
z  axis  is  at  the  starting  point  of  the  new  line  segment.  Similarly,  Eq.  71  translates 
the  time  axis  such  that  the  starting  time  of  the  polyline’s  new  segment  is  set  equal 
to  the  ending  time  of  the  first  segment,  figure  5  illustrates  the  result  of  applying 
Eqs.  70  and  71  to  Eig.  3a. 

Our  two  RES  polyline  algorithms,  R1.S1  and  RI  S2,  are  distinguished  by  the  means 
used  for  detecting  innovation  bias.  In  RLS1,  we  first  divide  y(k)  by  the  square  root 
of  its  covariance,  given  by  Eq.  63,  to  obtain  the  normalized  innovation,  S(k): 


S{kt)  +  j)  —  >(A„  +  j) 


I  /  ./  -  1  2j  -  I 

a  \  j  +1  2j  +  1 


Bias  detection  is  then  accomplished  by  applying  a  fixed  threshold,  t,  to  S(A'): 

If  | S(0  +  1)|  >  < 

(73) 

Then  terminate  current  segment  at  k  =  0  . 

As  we  will  discuss  later  in  Section  6,  selecting  a  threshold  to  provide  cither  a  fixed 
error  norm  or  a  specified  data  compaction  ratio  requires  iteration  on  the  value  oft. 

Our  bias  detector  RES2  is  composed  of  two  independent  criteria  such  that  a  bias 
is  declared  if  either  criterion  is  satisfied.  One  of  the  two  bias  detectors  in  RLS2  is 
fixed  threshold  detection  applied  to  ,S’(A),  i.e.,  Eq.  73.  The  second  bias  detector  in 
RI.S2  is  as  follows: 

If  ■)(!)  -t  n)  ■  y(()  +  n  r  1)  >  0  ,  n  1,2,  ....  (N  -  I) 

(74) 

Then  terminate  current  segment  at  A  =  (I  . 


Figure  5  Application  of  Eqs.  70  and  71  to  Fig.  3a.  The  polyline  shown  in  Fig. 
3a  has  been  translated  such  that  the  end  of  the  first  line  segment  is  at  the  start¬ 
ing  point  of  the  second  line  segment. 
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In  words,  bias  is  detected  as  a  sequence  of  N  values  of  -y(^)  all  having  the  same 
sign,  where  /V  (like  e  in  Eq.  73)  is  a  threshold  parameter.  The  idea  behind  Eq.  74 
is  that  the  innovation  sequence  should  be  spectrally  white,  and  therefore  subject  to 
frequent  sign  changes,  only  so  long  as  the  one-segment  waveform  model  (Eq.  52) 
is  valid.  A  test  similar  in  spirit  to  Eq.  74  was  discussed  by  Pavlidis  (Ref.  7,  p.  288, 
Fig.  12.5). 

Our  discussion  of  the  performance  of  the  RLS1  and  RLS2  polyline  algorithms 
is  deferred  until  Section  6,  where  their  performances  are  compared  to  those  of  a 
number  of  other  methods. 


4.3  REFORMULATED  STATE  EQUATIONS 

In  this  section  we  slightly  recast  our  one-parameter  regression  results  into  a  form 
suitable  for  use  with  the  general  GLR  formalism  presented  in  Section  2. 

From  Eq.  46, 


P(k)  =  cov[jc(At)1  =  E  r  |  •  (,ys| 


_  fcov(X 
(  Ely  • 


|y(y  )l  Ely  •  s] 
Lv  •  s]  cov[s(/)] 


However,  from  Eq.  51a, 


cov[y]  =  E[y  •  s]  =  0  , 


and,  from  Eq.  64, 


cov[s(/j]  =  P,  (J)  = 


j  ■  (j+  1 )  ■  (  2j  +  I ) 


From  Eqs.  75  and  76, 


Pik)  =  P,U) 


n 


From  Eqs.  46  and  56,  we  simplify  notation  by  writing 


x(k„  +  j)  =  j^j  =  x(n\k„  + 


independent  of  n. 

From  Eqs.  12,  49,  77,  and  78, 

K(k)  =  Ps(j  -  1) 


n  (’] 


V  '(A) 


K{k)  =  [j  ■  P,(j  -  \)/V(k)\ 
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From  Eqs.  79  and  12 ' , 


}•; 


where  A,  (/)  is  given  h\  lie).  62. 


5.0  (>I.K  POLYLINK.  ALGORITHMS 


5.1  SI  MM\K\  (II  PRI  VIOl  S  RKSI  I  IS 


In  this  section  we  gather  together  tor  convenient  reference  the  various  previously 
derived  i|iKintiiics  needed  to  partieulaii/e  the  general  Cil  R  formalism  to  the  poly¬ 
line  approximation  problem. 

- (i  ijm 

■I1 1 A  t  I.A)  /  (47] 

. \  ') 


v'  V 

eg 
v  SS 


//(A,,  t  /)  i  I  /  | 

—  iV'Ti.1 


/’(A)  /’  (/!  ■  jj| 


x  <  A,,  r  /) 


V  ( »l  A,,  4  /) 


A  (A)  A  (/) 


(7  f  I )  -  (2/  4  I  ) 


1(A)  =  !,,(/)  a 


(!*!)(:;':) 
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j  ■  u  +  i)  •  (2j  +  i) 


We  recall  that,  throughout  the  above  equations. 


j  =  k  —  kt  |  , 
r  ss  0  -  kn  . 


From  Eqs.  16-18  and  47  we  have 


k  <  6 
k  >  0  ' 


5.2  THE  AUXILIARY  VARIABLE 


From  Eqs.  49  and  80, 


K(k)  ■  H(k)  =  /C,  O') 


[?  ?}• 


It  then  follows  from  Eqs.  81  and  82  that  the  difference  equation  for  the  auxiliary 
variable,  Eq.  25,  presently  simplifies  to 

e«„  +;;«)  =  {-k\j)  ii-,  m 


With  the  definition 


0(A„  +  +  r) 


f0u(y»  0i2  0») 

[e;,  (j\r)  o::  (./» j  ’ 


it  follows  from  Eq.  83  that 

0.1  (7)  =  e„(7  -  1)  (85a) 

0,:  0)  =  0i:  U  ~  I )  (85b) 

0,,  0)  =  -A',  O')  •  '),,  0  -  1)  +  (I  -  j  ■  A  (7)1  ■  02i  0  -  1)  (86a) 

022  O)  =  -A\0)  '  0.:  0  -  1)  +  II  -  J  ■  A\  0)1  ■  022  0  -  1)  •  (86b) 

where  for  conciseness  we  temporarily  suppress  the  second  index  (r)  in  the  various 
0„,„.  We  find  from  Eqs.  27  and  85  that 

0, ,(/;/-)  =  I  ,  (87a) 

0|,(j;r)  -  0  .  (87b) 

and  then  from  Eqs.  86  and  87  that 

0;i  (77)  =  -A,  0)  +  [I  -  ./  •  A,  0)1  ■  ();,(./  -  l;r)  .  (88a) 


K wyVr 


ss*; 


V.VJV.W.  ,  , 


w 


WJ 


WWW 
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e- (/';/•)  =  (I  -  j  ■  A.  (./)  1  ■  o r(y  -  I;/-)  .  (88b) 

The  solutions  of  Eqs.  88,  with  initial  conditions  given  by  Eq.  27,  are  developed 
in  Appendix  C  as 


0:,  (./»  = 

O;:  (j\r)  = 


3k(r+  1)  -  yq  +  1)1 
/(./  +  1)  (2j  +  1) 

r(/-  +  I)  (2/-  +  1 ) 
J(/V\)  (2j  +  1) 


(89a) 

(89b) 


Equations  87  and  89  provide  the  desired  closed  forms  for  the  elements  of  0. 


5.3  JUMP  SIGNATURES 

From  Eqs.  24,  81,  82,  and  89,  we  obtain 


phi  Us)  rr-Uy))  =  y  K  in)  .  [  1  0  ]  .  (0  0] 

(./»  /  (./;/•)]  hr  (0:,  6;:  j  U  n) 

=  £  A’,  («)  •  0::  (n\r)  •  . 


From  Eqs.  62  and  89b, 


6tt 


K  (n)  ■  ()„(/!;/•)  = 

J(j  4  1)  (27  +  1) 


From  Eqs.  90  and  91, 

Am  (y»  =  (./';/  )  =  0 

and 

6 


However,  since 


(90) 


(91) 


(92) 


A'„,  (>;/•)  =  -  -  •  V  n"‘  ,  //;  *=  1,2  . 

./  (y  +  1 )  (2j  +  1 )  h 


(93) 


and 


1 

,  ■  A/  +  1 )  . 


i 

■  ./(./  +  1)  (2/  +  1)  . 

6 


(94a) 


(94b) 


»  1 

x;  /(n)  =  x]  /(">  -  x;  /(«) . 

»  <1  «l 


(94c) 
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it  follows  from  Eqs.  92-94  that 
^ii  U',r)  =  0  , 

F i2  0;r)  =  0  , 


^21  U\r) 


Fii  (J;r) 


3 

j(j  +  1 )  (2 ]  +  1 ) 
1 

Ju  +  1)  (2 j  +  1) 


I/O  +  1)  -  r(r  -  1)  , 

I/O  +  1)  (2/  +  1)  -  r(r  -  1)  (2r  - 


(95a) 

(95b) 

(95c) 

1 )]  ■  (95d) 


From  Eqs.  23,  49,  and  81,  it  can  be  shown  that 

G,  ( j\r )  =  1  -  j  •  F,,  0  -  l;r)  , 
G2C/»  =  7  •  [1  -  F220  -  I;r)I  ■ 
It  follows  from  Eqs.  95  and  96  that 


G,  0» 


G2  C/;r) 


3r(r  —  1)  —  Q—  1)0+1) 
0  -  1)  (2y  -  1) 
r(r  -  1)  (2 r  -  1) 

U  ~  I)  (2y  -  1) 


(96a) 

(96b) 


(97a) 

(97b) 


Equations  95  and  97  are  the  desired  closed  forms  for  the  elements  of  Fand  C, 
applicable  when  j  2  r;  F  and  G  are  identically  zero  when  j  <  r  (cf.  Eq.  22). 


5.4  STATE  ESTIMATE  AND  INNOVATION  JUMP  RESPONSES 

As  noted  in  connection  with  Eqs.  20  and  21,  the  state  estimate  and  innovation 
developed  by  the  Kalman  filter  can  be  decomposed  as 

x(k)  =  x,(Jc)  +  .v,  (y';r)  ,  (98a) 

y(k)  =  yt(k)  +  y2  (j\r)  ,  (98b) 

where  we  have  defined 


*2  (J\r) 

=  F(k;6)  ■  v  , 

(99a) 

72  0';f) 

=  G(k:0)  ■  v  . 

(99b) 

We  refer  to  jE2  ar|d  72  as  the  state  estimate  jump  response  and  innovation  jump 
response,  respectively. 

We  obtain  from  Eqs.  99,  95,  97,  78,  and  48,  for  the  state  estimate  jump  response. 


■*2  0;r) 


a(J  -  r)  [(J  +  1)  -  r]  [(2 j  +!)+/•] 
JU  +  1 )  (2j  +  I ) 


(100) 
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and  for  the  innovation  jump  response, 


...  oc  ■  r(J-  -  r) 

"y->  (/;/")  = - - - ■ - ,  j  >  r 

U  -  1)  (2j  -  1) 


The  jump  responses  are  related  by  the  equation 
72  U\r)  =  z2  C /)  -  z2  (./» 

=  a  ■  (J  -  r)  -  j  ■  s2(J  ~  Ur)  • 
The  asysmptotic  behavior  of  y2  *s  derived  from  Eq.  101  as 


lim  7:  (/;/")  =  a  ■  r/2 


lim  7,  ( j;r )  =  a  ■  (J  -  r)  .  (104) 

>-» 

(2  r)/j-  0 

To  verify  the  correctness  of  Eq.  101,  we  operated  on  a  two-segment  polyline  (Fig. 
6a)  with  our  recursive  regression  algorithm,  Eqs.  60  and  61.  As  shown  in  Fig.  6b, 
the  difference  between  our  closed-form  expression  for  y2,  Eq.  101,  and  the  recur¬ 
sively  developed  innovation,  Eq.  61,  is  imperceptible  on  the  scale  of  the  figure. 
Equation  101  indicates  that,  in  general, 

sgn(72)  =  sgn(a)  ,  (105) 

'  i. .  Mgn  of  7:  is  the  same  as  the  sign  of  a.  In  Fig.  7  we  provide  an  example 
negative  a.  The  dashed  lines  on  Figs.  6b  and  7b  are  the  asymptotes  given  by 
qs.  103  and  104. 

Cit  u  numerical  experience,  McAulay  and  Denlinger1  have  assumed 

y2(J\r)  =  K  ■  (j  -  r)2  ,  (106) 

with  K  an  unspecified  constant,  independent  of  j  and  r.  (We  have  cast  McAulay 
and  Denlinger’s  Eq.  19  into  our  notation  to  facilitate  comparison  with  our  Eq.  101 .) 
We  have  not  explored  the  conditions  under  which  McAulay  and  Denlinger’s  ap¬ 
proximation  for  7;  will  provide  similar  results  to  our  exact  expression,  Eq.  101. 


5.5  INVERSE  COVARIANCE  OF  INNOVATION  ESTIMATE 

In  this  section  we  obtain  closed-form  expressions  for  the  elements  of  the  inverse 
covariance  of  P{k),  C(k\6),  originally  defined  in  Eq.  30. 

With  the  definitions 

C<*„  +M.+r>  -  gjg’}  (m 


G{k„  +  j\ka  +/•)»* 


=  fGi  U\r)'\ 
[G2U\r) j  ’ 


LV 


Ml 

Y,» 

5$ 

•**!» 


S  J: 


V" 


'v 
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o  *  C2{j'S)  =  (2r  -  1)  •  5,  -  r(r  -  1)  (2r  -  l)3  •  s2  , 
a2  ■  C}  (j;r)  =  i  •  (2r  -  l)2  ■  [s,  -  4r2(r  -  l)2  ■  s;]  , 
where  the  quantities  y  are  given  by 


s,  O'if)  -  r2(r  -  l)2  ■  £  l/(n2  -  1) 

n~r 

1  ,  1W,  ,,  r2(r  —  1 ) 2 ( 2/  -t-  1) 

=  -  •  r(r  -  1)  (2r  -  1)  -  - — — - — — — 

2  2j(j  +  1) 


(110b) 


SiU\r)  m  £  1/ (4«2  -  1) 

n  =  r 

_/ _ r  -  1 

2j  -  1  2r  —  1  ’ 

s  S  «2/(4w2  -  1) 


(111b) 


U  -  r  +  1 )  +5,  (/»  ■ 


Equations  1 10  and  1 1 1  provide  the  desired  closed  forms  for  the  elements  of  C(k\0). 
Although  C(k;6)  can,  in  principle,  be  computed  off  line  and  stored,  it  is  actually 
not  necessary  to  do  this,  as  we  will  show  in  Section  5.7. 


5.6  THE  FILTERED  INNOVATION 

In  this  section  we  provide  simple  recursive  filters  for  developing  the  elements  of 
d(k\8),  the  filtered  innovation,  originally  defined  in  Eq.  29. 

With  the  definition 


d(k0  +  j;k0  +  r) 


we  find  from  Eqs.  29,  63,  and  108  that 


=  (d\  US)! 
[d2U\r) j  ’ 


ct2  ■  d,  (y»  =  £  y(k„  +  n)  ■  [a2  •  G,  {n\r)/V0(n)] 


o1  ■  d2(j;r)  =  XI  7(*n  +  n)  ■  { a 2  ■  G2(n;r)/Vt)(n)] 


Substituting  Eqs.  97  for  G,  and  G,  and  Eq.  63  for  K0  in  Eq.  113,  we  obtain 


a1  ■  d]  \j;r)  =  £  y(k0  +  n)  ■  A(n\r) 

n  =  r 
) 

a 2  ■  d2(j's)  =  £  7(*<>  +  »)  '  B(n\r) 
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where 


3 r(r  -  1)  -  (/r  -  1) 

A(n;r)  =  - 

(//  +  1)  (2/7  +  1) 

,  r(r  ~  1)  (2r  -  1) 

B(n\r)  as - 

(n  +  1)  (2/j  +  1) 

With  the  definitions 

A  U‘J)  =  o2  ■  di  O', 7  -  i  +  1) 

D2  (/';/)  =  O2  ■  d2  0,7  -  i  +  1) 
it  follows  from  Eq.  114  that 


(115) 


U16) 


A  (j;i)  =  D,  0  ~  1;'  -  1)  +  A(j;j  -  i  +  1)  •  y(k0  +  j) 


D2  (J',i)  =  02(J  -!;/-!)  +  SO;./  -<+!)•  y(k„  +  j) 


017) 


Equations  117  are  the  desired  recursion  relations  for  generating  the  filtered  inno¬ 
vation  sequences  dt  and  d2.  Equations  1 17  arc  used  as  follows.  When  a  new  mea¬ 
surement  is  received  at  time  j,  the  one-parameter  regression  generates  via  Eq.  61 
a  new  innovation  value,  y(k„  +  J).  Our  interest  is  in  obtaining  the  values  of  dt  (j;r) 
and  d2{j\r)  for  the  M  “most  recent’’  values  of  r,  where,  to  make  our  example  con¬ 
crete,  we  select  M  =  20  (cf.  Fig.  2b).  Thus,  our  interest  is  in  obtaining  tf,  (j\r)  and 
d2(j;r)  for  r  =  j,  (J  -  1),  ....  (/  -  19).  With  Eq.  116,  we  formulate  our  problem  as 

Given  (AC/  -  1;/)  ,  i  =  1,2 . 20)  and  y(k„  +  j) 

018) 

Calculate  |A  f /;/)  ,  i  =  1,2,  ....  20) 

For  i  =  1,  we  obtain  from  Eq.  117 


A  O';  1)  =  A  (j  -  1;0)  +  A 07)  •  7(*<>  +  j)  ■  019) 


From  Eq.  1 16, 

A  0  -  1;0)  -  a2  ■  dt  (J  -  1 7)  ,  (120) 

and  from  Eqs.  22  and  29 

d]  {j-r)  =  0  ,  j  <  r  , 

so  that 

A  0  -  1:0)  =  0  .  (121) 
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The  solution  to  our  problem,  Eq.  118,  follows  from  Eqs.  117,  119,  and  121: 
A  0;1)  =  A(j'J)  •  7 Ao  +  j) 

A  C/;2i  =  <4  0;v  -  i)  ■  7(*„  +  j)  +  AC/  -  i;D 

(122) 

A  A3)  =  A  (A/  -  2)  •  7(^0  +  j)  +  At/'  -  1;2) 


A  O';20)  =  j4C/V  -  19)  •  7A0  +  j)  +  A  O'  -  1:19) 

The  solution  for  A  is  obtained  by  making  the  replacements 

A  -  A 

B  -  A 

in  Eq.  122. 

Although  the  coefficients  A(n\r )  and  B(n\r)  given  by  Eq.  115  could  in  principle 
be  calculated  off  line,  we  chose  in  our  implementation  of  the  GLR  algorithm  to 
alleviate  storage  requirements  (at  some  cost  in  execution  speed)  by  calculating  the 
coefficients  as 

A(n\r)  =  Ra(r)  ■  Nh(n )  -  N„(n) 

(123) 

B(n\r)  =  Rh  ( r )  ■  Nb{n) 

The  vectors  Ra ,  Na,  Rh,  and  Nh  were  calculated  off  line  and  made  available  to 
the  algorithm,  where 

/UO-  3r(r  -  1 )  ,  (1 24a) 


2  n  +  1 


(124b) 


Rh  (r)  =  r(r  -  1 )  (2r  -  1 )  , 


Nh(n)  = 


(n  +  1)  (2 n  +  1) 


( 1 24d) 


In  our  GLR  implementation  we  calculated  200  values  of  each  of  the  four  vectors 
defined  in  Eq.  124,  for  a  total  of  800  stored  coefficients.  Calculating  the  complete 
set  of  coefficients  A(rt'r)  and  B(n;r)  off  line,  rather  than  computing  them  as  needed 
via  Eq.  123,  would  have  required  storing  2  x  200  x  200  =  80,000  coefficients. 


5.7  GENERALIZED  LIKELIHOOD  RATIO 


From  Eqs.  37b,  50,  and  112, 


b(k0  +  j;k„  +  r)  =  d2  (j'J)  -  r  ■  dt  (J;r) 


.  •-  ■%.  r*  L*.  %  A  ^ 
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From  Eqs.  37a,  50,  and  107, 

a(k0  +  j;kt>  +  r)  =  r  •  Ct(j;r)  -  2 r  ■  C2(j\r)  +  C,  (y';r )  ■  (126) 

The  likelihood  ratio  is  then  obtained  from  Eq.  36  as 

l{k;6)  =  6:(A;0)/tf(A;0)  .  (36) 

At  every  time  increment  (i.e.,  when  k  increases  by  unity),  Eq.  36  is  used  to  com¬ 
pute  M  =  20  new  values  of  ft 

I \k;k  -  i  +  1)  ,  /  -  1,2,  ....  20  . 

Recursive  relations,  Eqs.  117,  are  used  to  establish  d,,  which  in  turn  enter  into  the 
numerator  of  Eq.  36  by  means  of  Eq.  125.  The  C,  coefficients,  Eqs.  1 10,  are  used 
in  off-line  calculations  of  a(k;6)  via  Eq.  126,  which  in  turn  enter  into  the  denomina¬ 
tor  of  Eq.  36. 

The  GLR  algorithm  flow  of  operations  is  shown  in  Fig.  2b. 


6.0  NUMERICAL  EXPERIENCE 


6.1  INTRODUCTION 

In  this  section  we  present  some  numerical  results  illustrating  the  effectiveness  of 
several  polyline  algorithms.  The  algorithms  evaluated  numerically  are  listed  in  Ta¬ 
ble  2,  grouped  according  to  speed  of  execution. 

Performance  in  all  cases  was  based  on  approximation  of  the  infrared  cloud/sky 
waveform  shown  in  Fig.  8a.""  The  various  regions  of  the  waveform  (broken 
clouds,  blue  sky,  etc.)  were  identified  by  eye  and  labeled  manually.  However,  one 
objective  of  the  current  work  is  to  facilitate  the  development  of  algorithms  capable 
of  automatically  segmenting  this  type  of  waveform  data.  An  illustrative  polyline  rep¬ 
resentation  of  the  data  is  shown  in  Fig.  8b,  as  provided  by  algorithm  GL.R-M  (dis¬ 
cussed  in  Section  6.5). 

Comparison  between  generally  similar  waveforms  (c.g..  Figs.  8a  and  8b)  is  facili¬ 
tated  by  the  use  of  numerical  measures  of  similarity,  or  “error  metrics.”  Our  selec¬ 
tion  of  error  metrics  (for  approximations  with  a  fixed  number  of  knots,  i.e.,  for 
fixed  data  compaction)  includes: 


;:K.  A.  Steinberg  and  M.  .1.  McHugh,  An  Error  Detection  and  Smoothing  Algorithm  for  Infrared 
Data,  IHU/API  TO  1355  (Apr  1986). 

;,1  .  M.  Howscr,  H  ide  Area  Guidance  and  Control  Program:  Investigation  of  Scanning  IK  Seek¬ 
er  Performance  in  Hack  ground  Clutter ,  IHU/API.  TO  1360  (Dec  1986). 


structure 
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Table  2 

Polyline  algorithms  evaluated  numerically. 


Algorithm 

Author/ Reference 

Method 

EMK 

McLaughlin  and 

Zacharski" 

E-maximal  knots 

DP 

Dunham"  and 

Papakonstantinou 1,1 

Dynamic  programming 

Gl.R 

Steinberg  (Sect.  5) 

Generalized  likelihood  ratio 

LSBPT 

Tomek: 

Longest  segment  between 
parallel  tangents 

CIM 

Williams4  and 

Sklansky  and  Gonzalez5 

Cone  intersection  method 

HOP-F 

Pavlidis7 

Hop-along  algorithm  (fast) 

HOPS 

Pavlidis7 

Hop-along  algorithm  (slow) 

BAD 

Wall  and  Daniellson 

Bounded  area  deviation 

RLS1 

Steinberg  (Sect.  4,  Eq.  73) 

Recursive  least  squares  (fast) 

RLS2 

Steinberg  (Sect.  4,  Eq.  74) 

Recursive  least  squares  (slow) 

6.2  A  NOTH  ON  KRROR  METRICS  AND  PERFORMANCE  COMPARISONS 

When  the  infrared  camera  used  to  obtain  Tig.  8a  is  operated  under  closed-cover 
conditions  (with  the  lens  cap  on),  the  output  waveforms  are  reduced  to  zero-mean, 
unit-variance  Gaussian  noise.  In  fact,  all  data  obtained  with  this  camera,  including 
Fig.  8a,  have  an  additive  Gaussian  noise  component  of  one  count  root  mean  square. 
In  this  section  we  discuss  how  additive  Gaussian  noise  of  known  level  affects  the 
relative  usefulness  of  alternative  error  metrics. 

It  is  common  in  the  literature  to  select  peak  error,  £'«,  as  the  primary  measure 
of  polyline  algorithm  performance.  Defining  K  as  the  number  of  knots  in  the  ap¬ 
proximation  and  P  as  the  number  of  points  in  the  unapproximated  digital  wave¬ 
form,  we  obtain 


Jim  =  0  ,  (126) 

k  -n 


a 


>  A* 


§ 


assuming  that  knot  locations  are  required  to  coincide  with  the  waveform’s  sample 
values,  a  requirement  commonly  imposed  by  designers  of  algorithms  that  minimize 
Ex  .  At  fixed  data  compaction  (i.e.,  for  a  fixed  value  of  K),  we  can  use  £*.  to  com¬ 
pare  the  goodness  of  fit  provided  by  alternative  algorithms;  the  algorithms  are  ranked 
according  to  how  close  Ex  comes  to  the  zero  ideal  value. 

An  alternative  approach,  which  we  adopt  in  the  present  work,  is  to  select  £,  as 
the  primary  performance  metric.  Ideal  performance  is  expressed  as 

£,  I  ,  (127) 

because  the  data  are  known  to  contain  an  additive  component  of  zero-mean,  unit- 
variance  Gaussian  noise.  If,  corresponding  to  approximations  A  and  B,  we  have 
root  mean  square  fitting  errors  £\  ,  and  £,,,,  and  if 
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0  <  <  £,„  <  1  .  (128) 

we  recognize  the  B  approximation  as  being  better  than  the  A  approximation;  both 
approximations  are  under-smoothing  the  data  and,  very  likely,  using  more  knots 
than  are  really  needed  to  represent  the  underlying  structure. 


6.3  AN  EXPERIMENT  AT  7:1  DATA  COMPACTION 

The  ten  algorithms  listed  in  Table  2  were  used  to  develop  50-knot  approxima¬ 
tions  to  Fig.  8a,  corresponding  roughly  to  7;  1  data  compaction.  The  results  are  sum¬ 
marized  in  Table  3. 


Table  3 

Performance  of  polyline  algorithms  in  developing  50-knot  approximations  to  an 
infrared  cloud/sky  waveform,  Fig.  8a  (data  compaction  -  7:1). 


$ 

IS 


Error 


Execution  E«  £,  £,  E^ 


Algorithm 

Time  (s) 

Average 

Avg.  Abs. 

rms 

Peak 

EMK 

34.02 

-0.81 

2.95 

3.74 

8.80 

DP 

12.14 

-0.36 

2.95 

3.69 

8.10 

GLR 

1.79 

-0.05 

2.24 

3.00 

11.24 

LSBPT 

0.24 

-  1.68 

4.26 

4.87 

11.75 

CIM 

0.35 

-0.81 

2.95 

3.74 

8.80 

HOP-F 

0.24 

0.23 

2.61 

3.59 

11.62 

HOPS 

0.54 

0.13 

2.60 

3.57 

10.70 

BAD 

0.17 

-  1.82 

4.21 

4.95 

13.77 

RLS1 

0.14 

-0.01 

2.76 

3.40 

9.26 

RLS2 

0.21 

-0.01 

2.76 

3.40 

9.26 

si 
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Inspection  of  Table  3  shows  that  DP  had  the  smallest  peak  error,  while  GLR  had 
the  smallest  root  mean  square  and  average  absolute  errors.  RLS  and  GLR  had,  by 
far,  the  smallest  average  error.  The  two  fastest  methods  were  RLS1  and  BAD,  with 
RLS  enjoying,  relative  to  BAD,  a  small  advantage  in  speed  and  large  advantages 
in  every  measure  of  fitting  accuracy. 

Although  McLaughlin  and  Zacharski  claim  that  EMK  “...allows... one  to  optimally 
approximate  the  data...”,M  it  appears  from  Table  3  that  this  is  not  so.  Although 
EMK  yields  excellent  results,  its  approximation  is  generally  very  similar  to  that  provid¬ 
ed  by  CIM,  while  CIM  executes  about  100  times  faster. 

The  slow  and  fast  variants  of  RLS  were  found  in  this  case  to  yield  the  same  ap¬ 
proximation,  although  in  our  later  tests  this  was  found  not  to  be  generally  true. 

Tomek’s  algorithm,  LSBPT,  is  of  academic  and  historical  interest  as  perhaps  the 
first  linear-time  polyline  algorithm  derived  from  geometrical  reasoning  (as  contrast¬ 
ed  with,  for  example,  the  statistical  and  mathematical  frameworks  of  least-squares 
and  spline  approximations).  While  LSBPT  is  very  similar  in  structure  to  the  subse¬ 
quently  developed  CIM,  Table  3  indicates  that  LSBPT  executes  about  50%  faster 
than  CIM  but  provides  significantly  poorer  fitting  accuracy.  RLS  is  both  faster  and 
more  accurate  than  LSBPT. 

I  note  that  Tomek2  describes  both  a  fast  and  a  slow  variant  of  his  algorithm. 
Both  of  Tomek’s  algorithms  were  implemented,  and  it  was  found  that  the  slow  al- 
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gorithm  offered  almost  no  improvement  in  filling  accuracy.  I  consequently  report 
only  the  results  obtained  with  Tomek’s  fast  algorithm. 

The  fast  variant  of  Pavlidis’s  hop-along  algorithm,  HOP-F,  offers  significantly 
faster  speed  than  HOP-S,  for  a  relatively  small  penalty  in  fitting  accuracy. 

Our  preliminary  assessment  based  on  Table  3  is  as  follows: 

1 .  DP,  while  very  slow,  provides  smallest  peak  error,  E„ . 

2.  GLR,  provides  the  smallest  root  mean  square  error,  £,*  and  average  abso¬ 
lute  error,  £, . 

3.  RLS  provides  very  large  improvements  in  speed  relative  to  GLR  and  DP,  at 
relatively  low  cost  in  degraded  fitting  accuracy. 


6.4  A  CLOSER  LOOK  AT  THE  GLR  APPROXIMATION 

The  statistics  given  in  Table  3  are  useful  as  an  overall  performance  summary. 
In  this  section  we  present  additional  ways  of  examining  performance  intended  to 
provide  a  more  detailed  and  intuitively  meaningful  picture  of  algorithm  strengths 
and  weaknesses.  We  confine  the  discussion  in  this  section  to  GLR;  similar  details 
on  the  performance  of  several  other  algorithms  are  provided  in  Appendix  D. 

In  Fig.  9a  we  show  the  original  waveform  with  the  polyline  knots  superimposed 
(each  knot  is  denoted  with  a  +  symbol).  While  the  polyline  approximation  actually 
comprises  a  set  of  line  segments  connecting  the  knots,  we  have  elected  not  to  show 
the  line  segments  in  order  to  avoid  excess  detail  in  the  figure.  In  Fig.  9b  we  show 
the  pointwise  difference  between  the  original  waveform  and  the  polyline  approxi¬ 
mation,  i.e.,  the  pointwise  fitting  error.  The  histogram  of  the  pointwise  error  (Fig. 
9c)  is  nearly  Gaussian,  with  several  outliers  in  one  of  the  tails  of  the  density. 

Since  it  is  common  in  the  literature  to  see  emphasis  on  the  £„  metric,  we  present¬ 
ly  consider  the  performance  of  GLR  where  its  peak  error  is  worst,  viz.,  in  the  neigh¬ 
borhoods  identified  in  Fig.  9  by  the  circled  numbers  1-3.  These  peak-error 
neighborhoods  are  expanded  (zoomed)  in  Fig.  10;  Table  4  presents  the  data  from 
Fig.  10  in  tabular  form. 

Table  4 

Data  in  the  neighborhoods  of  the  three  points  of  worst  fit  for  the  GLR  algorithm 
(data  are  plotted  in  Fig.  10). 


s  5 

v,  JY 

» 

V,  $ 

^  s' 

V 

V 

'V  s 

I 


k 

z(k) 

Knot 

z(k) 

(z  -  z) 

Max.  Error  Vi 

22 

1336 

1335.91 

0.09 

23 

1296 

1298.60 

-  2.60 

iT-’ 

24 

1263 

8 

1261.30 

1.70 

25 

1254 

9 

1254.00 

0.00 

26 

1267 

1278.24 

-11.24 

^  f 

27 

1293 

1302.47 

-  9.47 

// 

28 

1325 

1326.71 

-  1.71 

u 

29 

1355 

1350.95 

4.05 

30 

1379 

10 

1375.18 

3.82 

275 

1106 

1108.48 

-  2.48 

276 

1115 

35 

1 1 1 1 .64 

3.36 

277 

1128 

1131.05 

-  3.05 

.  , 

278 

1152 

36 

1 1 50.47 

1.53 

279 

1200 

1211.18 

-  11.18 

-2 

280 

1269 

1271.89 

2.89 

y, 

40 

V 

THE  JOHNS  HOPKINS  UNIVERSITY 

APPLIEO  PHYSICS  LABORATORY 

LAUREL,  MARYLAND 


lIVIUUV. 


Table  4  (Continued) 

Data  in  the  neighborhoods  of  the  three  points  of  worst  fit  for  the  GLR  algorithm 
(data  are  plotted  in  Fig.  10). 
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6,5  A  HEURISTIC  MODIFICATION  TO  GLR 

Pernaps  the  most  important  observation  that  can  be  made  eoncerning  Fig.  9a  is 
that  GI.R  unnecessarily  uses  a  large  number  of  knots  to  represent  the  unstructured 
blue  sky  portion  of  the  waveform  (120  <  k  <  230).  By  representing  this  benign 
portion  of  the  data  more  efficiently,  we  can  free  up  additional  knots  needed  to  rep¬ 
resent  the  data  more  accurately  in  the  highly  structured  parts  of  the  waveform.  This 
suggested  that  we  develop  a  patch  to  our  GI.R  routine  (i.c. ,  additional  code  to  de¬ 
tect  the  blue  sky  region)  and  implement  the  RLS2  algorithm  rather  than  GLR  in 
this  region. 

We  observe  from  Fig.  9a  that  the  three  largest  errors  are  located  on  steep  shoul¬ 
ders,  i.e. ,  regions  of  the  waveform  where  the  slope  is  relatively  large.  This  suggests 
the  following  second  ad  hoc  modification  to  our  GLR  algorithm:  when  the  slope 
is  larger  than  some  threshold  value,  replace  the  nominal  likelihood  threshold  ( 
(Fig.  2)  by  a  smaller  number  a  ■  t.  where  0  <  k  <  1.  i.e.. 


.<•(./  )  <  f, 
V (./ )  >’(, 


where  s{j)  is  the  current  estimate  of  line  slope,  and  j  =  ( k  -  A„)  is  the  distance 
from  the  previous  knot.  An  appropriate  value  for  (,  may  be  developed  from  the 
data  in  Table  4. 

We  refer  to  the  resulting  modified  version  of  GLR,  which  incorporates  both  ad 
hoc  patches,  as  GLR-M.  The  50-knot  approximation  obtained  with  GLR-M  is  shown 
in  Fig.  11,  analogous  to  our  earlier  Fig.  9  results  for  GLR.  Comparing  Fig.  1 1  a 
with  Fig.  9a  we  note  that  the  blue  sky  region  is  represented  more  efficiently  by  GLR- 
M;  comparing  Fig.  1  lc  with  Fig.  9c  we  see  that  the  fitting  error  histogram  is  more 
compact  for  GI.R-M  than  for  GLR. 
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Figure  9  GLR  polyline  algorithm  evaluation,  (a)  IR  cloud  waveform  (solid  curve, 
z)  and  polyline  knots  ( + ).  (b)  Pointwise  difference  between  z  and  polyline  approx¬ 
imation.  (c)  Histogram  of  (b)  with  fitted  Gaussian  density  (circles  and  dashed  curve). 


In  Table  5  we  compare  GLR-M  ro  DP  and  lo  our  original  GLR  routine.  Com¬ 
pared  with  the  other  algorithms,  GLR-M  displays  significantly  superior  performance 
with  respect  to  £j  and  £; ,  while  simultaneously  achieving  a  value  of  Ex  just  4% 
worse  than  the  dynamic  programming  solution. 

The  patch  applied  to  GLR  was  tuned  specifically  to  achieve  a  50-knot  approxi¬ 
mation;  i.e.,  GLR-M  is  not  sufficiently  robust  to  be  directly  comparable  to  the  oth¬ 
er  algorithms  (which  is  why  GLR-M  results  were  not  included  in  Table  3). 
Nonetheless,  our  experiments  with  GLR-M  suggest  that,  with  additional  effort,  it 
may  be  possible  to  achieve  excellent  fitting  simultaneously  in  all  error  metrics. 
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Figure  10  Data  in  the  neighborhoods  of  the  three  points  of  worst  fit  for  the  GLR 
algorithm.  Figure  9a  shows  these  neighborhoods  within  their  larger  context;  Ta¬ 
ble  4  provides  these  data  in  tabular  form,  (a)  Neighborhood  of  peak  error  point, 
(b)  Neighborhood  of  second  worst  point,  (c)  Neighborhood  of  third  worst  point. 


Table  5 

Performance  of  three  polyline  algorithms.  Results  for  DP  and  GLR  as  for  Table 
3;  GLR-M  is  a  heuristically  modified  version  of  GLR. 
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Figure  11  GLR-M  polyline  algorithm  evaluation,  (a)  IR  cloud  waveform  (solid  cun/e) 
and  polyline  knots  (o).  (b)  Pointwise  fitting  error,  (c)  Histogram  of  (b)  with  fitted 
Gaussian  density  (circles  and  dashed  curve). 


6.6  ROBUSTNESS  WITH  RESPECT  TO  DATA  COMPACTION 

All  results  reported  in  Sections  6.3  through  6.5  have  been  for  50-knot  approxi¬ 
mations  to  Fig.  8a,  corresponding  to  a  fixed  =7:1  data  compaction.  In  this  section 
we  present  results  of  algorithm  performance  over  a  range  of  data  compactions,  viz., 
approximations  of  40-70  knots.  Our  approach  is  to  display  the  variation  of  each 
performance  metric  (En,  Et ,  £:,  £*, ,  and  T)  with  the  number  of  knots,  for  each 
algorithm. 
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Of  the  various  algorithms  described  m  the  litciaiuic,  CIM  provides  perhaps  the 
best  combination  of  good  fitting  accuracy,  fast  execution  speed,  robustness,  and 
ease  of  use.  For  these  reasons,  the  performances  of  the  other  algorithms  are  in  ev¬ 
ery  case  compared  with  that  of  the  CIM  baseline. 

Results  presented  in  this  section  compare  GLR  with  CIM  (Fig.  12)  and  RLS1  with 
CIM  (Fig.  13).  Analogous  figures  for  the  other  algorithms  are  provided  in  Appen¬ 
dix  E.  In  all  cases,  the  dashed  curves  are  for  CIM. 

The  plots  in  Figs.  12f  and  13f  refer  to  the  threshold  parameter  e  as  it  appears 
in  Eq.  73  for  RLS1  and  in  Fig.  2b  for  GLR.  The  parameter  €  for  CIM  is  simply 
the  peak  error: 


*  =  E„ 


for  GLR 


Thus,  as  previously  noted  by  Dunham,1  CIM  has  an  advantage  in  ease  of  use  rela¬ 
tive  to  algorithms  (such  as  RLS  and  GLR)  with  e  parameters  that  do  not  directly 
control  an  error  metric. 
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Figure  12  Performance  statistics  as  functions  of  data  compaction  for  GLR  (sol¬ 
id  curves)  and  CIM  (dashed  curves). 
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Figure  13  Performance  statistics  as  functions  of  data  compaction  for  RLS1  (solid 
curves)  and  CIM  (dashed  curves). 
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7.0  CONCLUSIONS 


The  extraction  of  shape  information  from  waveforms  can  be  facilitated  by 
preprocessing  the  waveform  with  a  polyline  algorithm,  i.e.,  an  algorithm  that  ap¬ 
proximates  the  waveform  as  a  concatenated  sequence  of  straight-line  segments.  Poly¬ 
line  approximation  accomplishes  smoothing  of  small-amplitude  structure  while  leaving 
large-amplitude  structure  well  defined  and  unsmoothed.  Moreover,  the  polyline  rep¬ 
resentation  generally  provides  a  large  degree  of  data  compaction  relative  to  the  original 
time  series. 

The  new  polyline  approximation  algorithms  developed  in  this  report  are  of  two 
types,  which  we  refer  to  as  recursive  least  squares  (RLS)  and  generalized  likelihood 
ratio  (GLR)  algorithms.  The  performance  of  these  algorithms  is  assessed  in  Section 
6  by  comparison  with  a  number  of  alternative  approaches.  The  Cone  Intersection 
Method  (C1M)  algorithm  discovered  independently  by  Williams4  and  by  Sklansky 
and  Gonzalez5  is  used  as  a  benchmark  against  which  to  compare  the  performance 
of  other  methods,  including  our  own  RLS  and  GLR  algorithms. 

The  relative  performance  of  the  various  algorithms  turns  out  to  be  somewhat  de¬ 
pendent  on  the  choice  of  error  metric.  In  particular,  the  original  algorithms  presented 
here  are  considerably  superior  to  prior  algorithms,  with  respect  to  bias  and  root 
mean  square  error. 

Our  RLS  algorithms  derived  in  Section  4  are  based  on  a  simple  Kalman  filter. 
Our  numerical  experience  is  that,  relative  to  the  best  of  the  prior  approaches,  RLS 
has  (Figs.  13  and  E-6) 

•  Faster  execution  speed 

•  Slightly  worse  peak  error 

•  Slightly  superior  root  mean  square  error 

RLS  enjoys  a  significant  speed  advantage  relative  to  the  other  algorithms  we  have 
evaluated. 

Our  GLR  polyline  algorithm  is  derived  in  Section  5  as  an  application  of  a  for¬ 
malism  originally  provided  by  Willsky  and  Jones. 1  Some  characteristics  of  GLR, 
relative  to  the  best  of  the  prior  approaches,  are  (Fig.  12) 

•  Worse  peak  absolute  error 

•  Superior  root  mean  square  error 

Although  most  of  our  results  indicate  a  trade-off  between  root  mean  square  fitting 
error  (£2)  and  peak  error  (£„),  our  numerical  experiments  with  a  currently  non- 
robust  GLR  variant  indicate  that  it  may  be  possible  to  obtain  excellent  fitting  in 
all  metrics  simultaneously. 

Conclusions  regarding  other  approaches  arc  as  follows: 

1.  The  method  of  £-maximal  knots'4  yields  nearly  identical  fitting  performance 
to  that  of  C1M  but  is  about  100  times  slower  in  execution  (Fig.  E-l  in  Appen¬ 
dix  E). 

2.  The  dynamic  programming  method  of  Refs.  10  and  11  provides  a  small  ad¬ 
vantage  in  £„  and  nearly  equal  £;  relative  to  CIM  but  is  about  30  times  slow¬ 
er  in  execution  (Fig.  E-2). 

3.  Tomek’s  algorithm2  executes  about  50%  faster  than  CIM  but  provides  sub¬ 
stantially  worse  fitting  accuracy  both  in  £,  and  in  f„  (Fig.  E-3). 

4.  Pavlidis’s  hop-along  algorithm  executes  about  50%  slower  than  CIM,  is 
slightly  superior  in  £,,  and  significantly  worse  in  £*  (Fig.  E-4). 

5.  The  bounded-area-deviation  algorithm  1,1  is  only  slightly  slower  than  RLS1; 
however,  it  provides  significantly  worse  values  of  £\  and  E,  than  RLS1  (Figs. 
E-5  and  E-7). 
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Appendix  A 

WILLSKY  AND  JONES  GLR  FORMALISM 


A.l  EXPLOITING  LINEARITY 

Willsky  and  Jones 17  provide  without  derivation  a  formalism  they  call  the  Gener¬ 
alized  Likelihood  Ratio  (GLR)  method  for  detecting  and  characterizing  discontinu¬ 
ous  (jump)  inputs  to  linear  systems.  The  elements  of  the  GLR  method  are  outlined 
in  Section  2.3  and,  for  completeness,  are  derived  here. 

We  assume  that  the  jump  detection  problem  is  formulated  as  per  Section  2.1. 
From  Eq.  1, 

x(k)  =  $(k,k  -  1)  ■  x(k  -  1)  +  T(A'  -  1)  •  w(k  —  1 )  +  8#k  ■  v  , 

(A- 1 ) 

from  which  it  is  obvious  that  x(k)  undergoes  a  jump  when  k  =  6. 

As  discussed  in  connection  with  Eqs.  15,  and  analogous  to  Eq.  98,  we  can  write 

x(k)  =  a,  (A)  4-  as  (A)  ,  (A-2) 

where  .v,  (k)  is  the  value  that  x(k)  would  take  in  the  absence  of  a  jump  (r  =  0  or 
k  <  0),  and  .y2(A)  is  the  state  perturbation  in  response  to  a  jump. 

As  first  noted  by  McAulay  and  Denlinger, 1  the  linearity  of  the  system  equations 
allows  their  solution  in  response  to  a  jump  excitation  to  be  written  as  the  sum  of 
a  jump-independent  stochastic  component  (a, )  and  a  second  component  (a2)  linear¬ 
ly  proportional  to  the  jump  amplitude.  Moreover,  with  the  jump  amplitude  and 
time  taken  as  deterministic  (though  unknown)  quantities,  the  jump-dependent  com¬ 
ponent  .v;  is  seen  to  be  deterministic. 

From  Eqs.  1,  2,  and  A-2, 


•v,  (A 

+ 

1) 

+ 

¥ 

II 

I, A )x,  (A) 

+ 

T (A)w(A) 

(A-3) 

zt  (A 

+ 

1) 

+ 

L< 

II 

1  )a'i  (A  + 

1) 

+  v( A  +  1) 

(A-4) 

and 

AS  ( A 

+ 

1) 

II 

& 

+ 

l,A)Xj  (A) 

+ 

5».*  +  i  ’  v 

(A-5) 

z2  (k 

+ 

1) 

=  H(k  + 

1  )as  ( A  + 

1) 

(A-6) 

Equation  A-5  is  subject 

to  the  initial  conditions 

AS  (d  - 

-  1)  =  o 

(A-7) 

A. 2  SOLUTION  FOR  a,  (A ) 


Setting  first  k  =  0  -  1  and  then  k  >  6  in  Eq.  A-5,  and  using  Eq.  A-7,  we  obtain 
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With  the  definition 


as  ( k )  =  $(*,0)  ■  v  , 


(A- 10) 


it  follows  from  Eqs.  A-7  to  A- 10  that 


*(*,0)  =  0  ,  k  <  0 


(A- 11) 


$(0,0)  =  / 

$(/(-  +  1,0)  =  $(£  +  1,*)  ■  $(A-,0)  , 

which  appeared  in  the  text  previously  as  Eqs.  16  to  18  (Section  2.3). 
Equations  A- 12  and  A- 13  can  be  solved  by  induction  to  obtain 


$(*  4-  1,0)  =  $(n  +  l.n) 


From  Eqs.  A-6  and  A-10, 


z2(k)  =  H(k)  ■  $(A:,0)  ■  p  . 


(A- 12) 
(A- 13) 


(A- 14) 


(A- 15) 


S: 


A. 3  JUMP  SIGNATURES 

The  linearity  of  the  Kalman  filter,  Eqs.  5-14,  allows  us  to  decompose  the  state 


estimate,  analogous  to  Eq.  A-2,  as 

x(k\n)  =  .y,  (Ar|«)  +  x2(k\n)  .  (A-16) 

From  Eqs.  7  to  9  and  13,  we  have 

x2{k  +  1 1 A' )  =  $(£  +  \,k)x2  ( A | A)  (A-17) 

x2(k\k)  =  as  (A|A  -  1)  +  K(k)y2(k)  (A-18) 

7 ,(A)  =  z2(k)  -  H(k)x2{k\k  -  1)  .  (A-19) 

From  Eqs.  A-I5  and  A-19, 


y2(k)  =  H(k)  ■  [*(*.0)1'  -  x2  ( A  |  A  -  D]  , 
which,  with  Eq.  A-17,  becomes 

y 2(k)  =  H(k)  ■  |*(A,0)i<  -  $(A-,A-  -  1  ).v2  ( A  —  1  |A  —  1 )]  .  (A-20) 
With  the  definitions 

,v:(A|A)  =  F(A;0)  •  i>  (A-21) 

and 
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Eq.  A-20  becomes 

G(k\8)  =  H(k)  ■  [*(*,0)  -  *(k,k  -  \)F(k  -  1;0)]  .  (A-23) 

Equation  A-23  appeared  previously  as  Eq.  23  in  the  text. 

It  can  be  shown  from  Eqs.  A-17,  A-18,  and  A-20  that 

*2  {k\k)  =  [I  -  K(k)H(k))$(k,k  -  \)x2(k  -  1  |Ar  -  1) 


+  K(k)H{k)*(k,d)v  . 


This  can  be  written  as 


F(kfi)  =  A(k)F(k  -  1;0)  +  B(k) 


where  we  define 


F(k;0)  •  v  a  x2  (k\k)  , 

A(k)  =  (/  -  K(k)H(k)]*(k,k  -  1) 
B(k)  =  K{k)H(k)*(k,6)  . 


(A -24) 


(A-25) 


(A-26) 

(A-27) 


(A-28) 


Defining  the  auxiliary  variable  0(/t;0)  as  the  homogeneous  solution  of  Eq.  A-25, 


Q(k\6)  =  A(k)0(k  -  1;0)  , 


0(0;0)  =  /  , 


(A-29) 


(A-30) 


it  can  be  shown  by  back  substitution  into  Eq.  A-24  that  the  solution  of  Eq.  A-25 
is  given  by 


F(k;d)  =  £  Q(k\n)B(n)  . 


(A-31) 


Equations  A-27  and  A-29  lead  directly  to  Eq.  25  in  the  text,  and  Eqs.  A-28  and 
A-31  lead  to  Eq.  24. 

From  Eqs.  A-24  and  A-26, 

F(k,Q)  =  K(k){H(k)[*(k,d)  -  <t>(k,k  -  1) F(k  -  1;0)]| 

+  $(k,k  -  1  )F(k  -  1;0)  , 
which,  with  Eq.  A-23,  may  be  written  as 


F{k;6)  =  K(k)G{k;6)  +  i>(k,k  -  1  )F(k  -  1  ;0) 
Equation  A-32  appears  as  Eq.  49  in  Ref.  17. 


(A-32) 


WWW 


V  V  V  '  *  "ji  *  *  "  . 
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A.4  THE  LIKELIHOOD  RATIO 

Jump  detection  is  formulated  as  a  hypothesis  testing  problem  as  follows:1 
H0 :  Hypothesis  that  no  jump  has  occurred,  i.e.,  that  k  <  6 

(A-33) 

H | :  Hypothesis  that  a  jump  has  occurred,  i.e.,  that  k  >  0  , 
and,  equivalently. 


Ha- 

y(k) 

=  >1  (At) 

Hr- 

y(k) 

=  Ti  ( k )  +  G(A';0) 

(A-34) 


where  -y,  (A)  is  zero  mean  and  spectrally  white. 

Maximum  likelihood  estimates  for  d  and  v,  denoted  6  and  v,  are  obtained  im¬ 
plicitly  as  the  values  of  6  and  v  that  maximize  the  joint  conditional  density: 


p["Y C 1 ) . y(k)\Hit6,v\  . 

The  choice  between  Ha  and  //,  is  based  on  the  criterion 

A(Ar)fi  I  Zh  >  (A'35) 

(.<  v  =*  "o 

where  t)  is  a  threshold  ultimately  relatable  to  probabilities  of  false  detection  and  missed 
detection,  and  the  generalized  likelihood  ratio  A {k,6)  is  defined  as 


Pi  7d) . 

Ph(  1).  y(k)\Hn] 


(A-36) 


Since  the  conditional  densities  in  Eq.  A-36  are  both  Gaussian,  we  can  simplify  the 
implementation  of  Eq.  A-35  by  taking  the  logarithm  of  both  sides,  to  obtain 

I’m  =  2  ln[A(M)l  =  !j  -  f.  ,  (A-37) 


where 

A 

<j  -  E  y'  ■  K  '  ■  7„  (A-38) 

n  -  w 

and 

A 

I;  -  E  -  G,l„  ■,■)’  ■  l'-1  ■  (7fl  -  G„,  ■  v)  .  (A-39) 

n-li 


For  conciseness  of  notation,  subscripts  and  arguments  are  used  interchangeably  in 
Eqs.  A-38,  A-39,  and  the  rest  of  this  section.  For  example,  GlH  =  G(j\0).  Also, 
we  recall  that  Vn  in  Eqs.  A-38  and  A-39  is  the  covariance  of  7  (cf.  Eq.  1 1  in  the 
text),  while  G„„  ■  v  in  Eq.  A-39  is  the  mean  value  of  7.  subject  to  hypothesis  //, . 


-v  ^  n.'  V  V  t  *-  V’  t- 
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Equation  A-39  can  be  expanded  to  obtain 


(>  =  £  yl  ■  K'  y„  +  £  vT  •  ( Gl ,  •  V?  ■  G„'t)  ■ 


-  £  hi  ■  Vn'  ■  G„  •  *  +  vT  ■  Gle  •  V-'  •  7„)  . 


From  Eqs.  A-37,  A-38,  and  A-40, 


(A-40) 


I'm  =  •  T  £  Git  ■  v;'  ■  G».»l  ■  »  +  f  £  Gle  •  [/;'  •  7  1  • 

+  " T  '  [  £  GI»  '  V«  1  •  T/i  j  • 


(A-41) 


With  the  definitions 


and  noting  that 


we  write  Eq.  A-41  as 


du  =  £  Glt  ■  V- 


Ck.e  *  £  Gl,  ■  V-'  ■  Gn.t 


dlt»  ■  v  —  vT  ■  dke  , 


=  -vT  •  Ck„  ■  v  +  2dl  ■  v  . 


(A-42) 


(A-43) 


(A-44) 


(A-45) 


Equations  A-42  and  A-43  were  given  previously  in  the  text  as  Eqs.  29  and  30, 
respectively. 

Maximum  likelihood  estimates  for  v  and  0,  denoted  v  and  0,  are  obtained  by  set¬ 
ting  equal  to  zero  the  appropriate  partial  derivative  of  Eq.  A-45,  as  follows.  We 
note  that 


V,(/  ■  CkB  ■  v)  =  2Cke  ■  v 


^ i. (dk,n  •  v)  —  dke  . 


From  Eqs.  A-45  to  A-47, 


V  „  (4,« !  =  -2Cm  •  «  +  2dk  „ 


(A-46) 


(A-47) 


(A-48) 
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However, 


v,(4») 


=  o 


From  Eqs.  A-48  and  A-49, 

*'*  =  C*.tf  •  dke  , 

which  appeared  in  the  text  as  Eq.  28. 

From  Eqs.  A-45  and  A-50, 

(u  =  -dT(C~')TC  C“'  d  +  2dT  C_l  d  . 


(A-49) 


(A-50) 


However,  since 


C  =  CT 


and 

C"‘  =  (C'l)T  , 


we  obtain 


f(*;fl)  =  dT(k\d)  ■  C'l(k;6)  ■  d(k;6)  ,  (A-51) 

which  appeared  in  the  text  as  Eq.  32. 


Appendix  B 

ONE-PARAMETER  RECURSIVE  REGRESSION 


In  this  Appendix  we  derive  our  one-parameter  recursive  regression  algorithm,  Eqs. 
8 '  — 14 ' . 

From  Eqs.  12'  and  14', 


p, 0)  =  [i  -r  •  ps’u  -  \)/v(k)]  ■  p,u  -  i) 

Substituting  Eq.  11'  into  Eq.  B-l,  where 


we  obtain 


V(k)  =  j  ■  PJj  -  \)  +  o'  , 


(B-l) 


(ID 


(B-2) 


58 


j; 

r>  v- 
0:1 
>■  '1 

» 

& 

£  AJ 

£ 

4 
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y(k)  =  z(k)  -  j  ■  s(J  -  1)  ,  (B-14) 

which  appeared  in  the  text  as  Eq.  61. 

Equation  60  requires  no  derivation,  having  appeared  previously  as  Eq.  13': 

s(j)  =  W  ~  1 )  +  KsU)  ■  y(k)  .  (B- 15) 

Equations  B-15,  B-14,  B-13,  B-12,  and  B-10,  appearing  in  the  text  as  Eqs.  60-64, 
provide  the  desired  one-parameter  recursive  regression  procedure. 


Appendix  C 

SOLUTION  FOR  THE  AUXILIARY  VARIABLE 


In  this  section  we  derive  Eqs.  89  as  the  solutions  of  the  difference  Eqs.  88  subject 
to  initial  conditions  given  by  Eq.  27. 

With  the  definitions 


(C-l) 

B, 

=  ©22  (/» 

(C-2) 

Cj 

=  K'U)  =  6/lU  +  l)(2j  +  1)1  , 

(C-3) 

Eqs.  88  are  written  as 


= 

0-7 

■  Cj)  ■ 

-  C, 

(C-4) 

B,  = 

0-7 

■  ct)  ■ 

S,-i 

(C-5) 

From  Eq.  C-3, 

(1  -  J 

f  • > 

u  - 

1 ) (2 j  -  1) 

(C-6) 

u  + 

1 )  (2y  +  1) 

From  Eqs.  C-5  and  C-6, 

U  +  D(2y  +  1) 

•  B,  = 

J(2j  - 

1)  ■  Br  ,  -  (2 j  -  1)  •  Bl  X 

.  (C-7) 

With  the  definition 

D,  - 

7  (2/  - 

1  )•»,,. 

(C-8) 

Eq.  C-7  is  written  as 

j  ■  di4 

i  =  (./ 

-D-D,  , 

60 


mm* 
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©21  (r-r)  =  0  , 


(C- 


Eq.  C-19  becomes 


©21 


3 \r(r  +  1)  -j{J  +  D) 

JU  +  1)  (2 j  +  1) 


which  appeared  in  the  text  as  Eq.  89a. 
which  has  the  solution 


(C- 


(y  —  1 )  •  D,  =  constant  .  (C 

From  Eqs.  C-2,  C-8,  and  C-9, 


©22  0';D 


constant 

y(7  +  I)(2y  +  1) 


Obtaining  the  appropriate  initial  condition  from  Eq.  27, 

©22  (r’r)  =  1  » 


(C- 


(C- 


Eq.  C-IO  becomes 


©22<j;r) 


r(r  +  1 ) (2 r  +  1) 

jrr+T)(2T+T)  ' 


(C- 


which  appeared  in  the  text  as  Eq.  89b. 

From  Eqs.  C-3,  C-4,  and  C-6, 

(7  +  1 )  ( 2/  +  1 )  ■  /4,  —  (7  —  0(27  -  1)  ■  Aj_  i  -  6  .  (C- 


Defining 


E,  =  (7  -  1 )  (2y  -  1 }  •  ,  (C- 


Eq.  C-13  becomes 


Defining 


(7  +  1)  •  i  =  j  ■  Ej  -  67  . 

/■;  -  j  ■  £,  . 


Eq.  C-15  becomes 


£,•  1  =  r, 


6/  ■ 


(C 


(C 


(C 


We  can  show  by  back  substitution  that  the  solution  of  Eq.  C-17  is  given  by 


>Vik>!iAM!kria* 
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F,  =  a  -  3/(y  -  1 )  , 


where  a  is  a  constant. 

From  Eqs.  C-l,  C-14,  C-16,  and  C-18, 


a  -  3/(7  +  1) 

9  (j-r)  =  - -  — - 

7(7  +  0(2/  +  ! ) 

Obtaining  the  initial  condition  from  Eq.  27, 


(C-l  9) 


Appendix  D 

NUMERICAL  EXPERIMENTS  AT  7:1  DATA  COMPACTION 


In  Section  6.4  we  presented  graphs  depicting  polyline  fit,  fitting  error,  and 
histogrammed  fitting  error,  for  the  GLR  polyline  algorithm.  Here  we  provide  results 
in  an  identical  format  for  the  following  additional  algorithms: 


EMK 

CIM 

DP 

LSBPT 

HOPS 

BAD 

RLS1 

RLS2 


50- knot 
50-knot 

49- knot 

50- knot 

51 - knot 
50-knot 
50-knot 
50-knot 


approximation 

approximation 

approximation 

approximation 

approximation 

approximation 

approximation 

approximation 


(Fig.  D-l) 
(Fig.  D-l) 
(Fig.  D-2) 
(Fig.  D-3) 
(Fig.  D-4) 
(Fig.  D-5) 
(Fig.  D-6) 
(Fig.  D-6) 


Not  all  algorithms  were  able  to  provide  exactly  a  50-knot  approximation;  conse¬ 
quently,  the  actual  numbers  of  knots  are  given  above  for  the  fit  shown  in  the  fig¬ 
ure.  This  also  explains  the  discrepancy  between  some  of  the  performance  statistics 
appearing  below  the  histograms  in  the  figures  and  the  corresponding  results  in  Ta¬ 
ble  3  (in  the  main  text).  If  a  particular  algorithm  could  not  obtain  the  desired  50-knot 
approximation,  the  statistics  in  Table  3  are  an  average  for  fits  of  49  and  51  knots. 

Fits  provided  by  algorithms  EMK  and  CIM  were  identical  in  this  instance  and 
so  are  provided  as  a  single  figure.  Similarly,  results  for  RLS1  and  RLS2  were  identi¬ 
cal  for  the  50-knot  case  and  are  therefore  represented  by  a  single  figure. 

References  to  the  literature  and  definitions  of  the  algorithm  acronyms  are  provided 
in  Table  2  (in  the  main  text). 
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Azimuth 


Bin  size  1.090  Total  counts  371 

Minimum  -8.193  Average  -0.314 

Maximum  8.161  Sigma  3.751 

Figure  D-2  DP  polyline  algorithm  evaluation,  (a)  IR  cloud  waveform  (solid  curve, 
z)  and  polyline  knots  (o).  (b)  Pointwise  difference  between  z  and  polyline  approxi¬ 
mation.  (c)  Histogram  of  (b)  with  fitted  Gaussian  density  (circles  and  dashed  curve). 
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1  3  5  7  9  11  13  15 

8 in  si^e  1.494  Total  counts  371 

Minimum  10.667  Average  1.678 

Maximum  11.750  Sigma  4.875 

Figure  D-3  LSBPT  polyline  algorithm  evaluation,  (a)  IR  cloud  waveform  (solid 
curve,  z)  and  polyline  knots  (o).  (b)  Pointwise  difference  between  z  and  polyline 
approximation,  (c)  Histogram  of  (b)  with  fitted  Gaussian  density  (circles  and  dashed 
curve). 
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Bin  si/v  1  360  To1.il  counts  371 

Minimum  10  400  Avorcicjo  0  164 

Maximum  10.000  5>ic|ma  3  397 


Figure  D-4  HOPS  polyline  algorithm  evaluation,  (a)  IR  cloud  waveform  (solid 
curve,  z)  and  polyline  knots  (o).  (b)  Pointwise  difference  between  z  and  polyline 
approximation,  (c)  Histogram  of  (b)  with  fitted  Gaussian  density  (circles  and  dashed 
curve). 


L*j 


*■#/»*,* 


I 

I 


THE  JOHNS  HOPKINS  UNIVERSITY 

APPLIED  PHYSICS  LABORATORY 

LAUREL.  MARYLAND 


1 

B 

E 


B 

I 

I 

I 

e 

B 
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Azimuth 


Bin  size  1.701  Total  counts  371 

Minimum  13.769  Average  1.825 

Maximum  11.750  Sigma  4.958 

Figure  D-5  BAD  polyline  algorithm  evaluation,  (a)  IR  cloud  waveform  (solid  curve, 
z)  and  polyline  knots  (o).  (b)  Pointwise  difference  between  z  and  polyline  approxi¬ 
mation.  (c)  Histogram  of  (b)  with  fitted  Gaussian  density  (circles  and  dashed  curve). 
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Bin  size  1.191  Total  counts  371 

Minimum  -8.604  Average  -0.010 

Maximum  9.256  Sigma  3.405 

Figure  D-6  RLS1  and  RLS2  polyline  algorithm  evaluation,  (a)  IR  cloud  waveform 
(solid  curve,  z)  and  polyline  knots  (o).  (b)  Pointwise  difference  between  z  and  poly¬ 
line  approximation,  (c)  Histogram  of  (b)  with  fitted  Gaussian  density  (circles  and 
dashed  curve). 
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Appendix  E 

PERFORMANCE  STATISTICS  VERSUS  DATA  COMPACTION 


In  Section  6.6  we  presented  curves  of  performance  statistics  versus  data  compac¬ 
tion  for  the  GLR,  RLS1,  and  CIM  polyline  algorithms.  Here  we  provide  results 
in  an  identical  format  for  an  additional  set  of  algorithms: 

EMK  (solid  curves)  versus  CIM  (dashed  curves)  (Fig.  E-l) 

DP  (solid  curves)  versus  CIM  (dashed  curves)  (Fig.  E-2) 

LSBPT  (solid  curves)  versus  CIM  (dashed  curves)  (Fig.  E-3) 

HOP-S  (solid  curves)  versus  CIM  (dashed  curves)  (Fig.  E-4) 

BAD  (solid  curves)  versus  CIM  (dashed  curves)  (Fig.  E-5) 

RLS2  (solid  curves)  versus  CIM  (dashed  curves)  (Fig.  E-6) 

RLSI  (solid  curves)  versus  RLS2  (dashed  curves)  (Fig.  E-7) 

References  to  the  literature  and  definitions  of  the  algorithm  acronyms  are  provided 
in  Table  2  (in  the  main  text). 

Of  the  various  algorithms  described  in  the  literature,  CIM  provides  perhaps  the 
best  combination  of  good  fitting  accuracy,  fast  execution  speed,  robustness,  and 
ease  of  use.  For  this  reason,  the  performance  of  the  other  algorithms  is  generally 
compared  with  that  of  the  CIM  baseline.  The  exception  to  this  is  Fig.  E-7,  where 
we  compare  the  two  recursive  least  squares  (RLS)  algorithms  developed  from  ana¬ 
lytical  considerations  in  Section  4. 
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Figure  E-3  Performance  statistics  as  functions  of  data  compaction  for  LSBPT 
(solid  curves)  and  CIM  (dashed  curves). 
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Figure  E-4  Performance  statistics  as  functions  of  data  compaction  for  HOPS 
(solid  curves)  and  CIM  (dashed  curves). 
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Figure  E-7  Performance  statistics  as  functions  of  data  compaction  for  RLS1  (solid 
curves)  and  CIM  (dashed  curves). 
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